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1 6. Abstract 


The overall objective of this study was to develop the analytical approaches that will lead to the solution \ 
of two classical image processing problems: 1) Image correlation and 2) Image sampling. Relevant multi- ' 

spectral image statistics which are applicable to image correlation and sampling are identified and developed. 

Task' I Image Statistics - In this task, relevant multispectral general image statistics and difference 
(change detected) image statistics were determined. The statistics of time sequential imagery in three 
spectral bands and representative of four general terrain types are presented. The general image statistics 
selected include intensity mean, variance, amplitude histogram, power spectral density function, and auto- 
correlation function. Difference image statistics presented include the intensity mean, variance, and 
quantized probability distribution, i 

Task II Image Correlation - This study task considered the translation problem associated with digital 
image registration and developed the analytical means for comparing commonly used correlation techniques. 

Using suitably defined constraints, an optimum and five suboptimum image registration techniques are 
defined and evaluated. A computational comparison was made. Both Gaussian and real difference image 
statistics were used to compare the techniques in terms of radial position location error. 

Task III Image Sampling - In this task, general expressions for determining the reconstruction error for 
specific image sampling strategies are developed. The image sampling parameters considered include 
scanning apertures, scanning lattices, image statistics, interpolation functions, and image recorder 
characteristics. These analytical expressions were programmed on a digital computer and the effects 
of sampling parameter variation demonstrated. ' ! 
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PREFACE 


The objective for this study was to develop the analytical approaches that 
lead to the solution of two classical image processing problems, image corre- 
lation and image sampling. Relevant multispectral image statistics applicable 
to image correlation and sampling were identified and developed. 

The study was divided into three separate but interrelated tasks. In the 
first task, relevant general image statistics and difference image (change 
detected) statistics were determined by suitable photographic processing 
and digital processing of multispectral imagery. The statistics of this 
time sequential imagery in three spectral bands (green, red, IR) and represen- 
tative of four terrain types (farmland, urban, mountains, forest) were presen- 
ted as numerical and graphical data. The second task considered the trans- 
lation problem associated with digital image registration and developed 
analytical means for comparing an optimum and five commonly used correlation 
techniques. A digital simulation was developed to compare the techniques in 
terms of radial position location error using both Gaussian and real difference 
image statistics. In the third task, general expressions were developed for 
determining digital image reconstruction error for specific image sampling '• 
strategies and sensor systems. Analytical expressions were programmed on a 
digital computer and the effects on reconstruction error caused by image 
sampling parameter variation were shown. 

Relevant image statistics for image correlation include the image mean, 
variance, autocorrelation function, and difference image histogram. Of 
these, the difference image histogram was found to be the most useful. From 
two statistical runs the tentative conclusions are: 

The optimum image registration algorithm appears to be a valid means for 
comparing the position location performance of commonly used techniques. 

On a computational basis, the mean absolute difference algorithm (using 
adaptive thresholding) and the cross correlation coefficient (using FFT 
techniques) are the most attractive correlation algorithms. 
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In general, image registration using hard limited (2 level) images 
resulted in slightly greater radial error than when using 6U level 
images . 

Those image statistics found to he relevant for image sampling are the 
power spectral density functions. The analytical approach was developed 
■to describe the effects of parametric variation of sample density, terrain 
type, recording device characteristics and presample filtering in terms of 
rms reconstruction errors. It was found that the mean squared error criterion 
does not restrict the choice of sampling aperture shape or size as a function 
of sampling lattice unless the sampling is very anisotropic. Computer image 
interpolation to a refinement of five recording device spots (in each direction) 
per pair of sensor sample values is sufficient. 

Future image correlation study should be directed toward extending the digital 
image correlation simulation. For image sampling, sensor system noise effects 
and means for reducing such effects merit attention. A promising area for 
investigation is the development of improved image reconstruction methods. 
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SECTION 1 

MULTISPECTRAL IMAGE STATISTICS 


1 . 1 INTRODUCTION 


In this section, multispectral imagery is selected and digitized to provide 
representative multispectral image data seen by an advanced spaceborne earth 
observation platform. General multispectral image statistics are determined 
which include image intensity mean, variance, amplitude histogram, power 
spectral density function array and autocorrelation function array. In 
addition, relevant difference or "noise" image statistics are generated 
which describe the "changes" that occur in imagery taken on sequential 
overflights. These statistics include the difference image intensity mean, 
variance, and quantized probability histogram. 

The analytical use of these statistics for performing image correlation and 
sampling studies is discussed in Section 2 and Section 3 respectively. 

1.2 IMAGE SELECTION 

In recent years, much progress has been made in the observation of major 
geographical areas in the continental United States using a variety of multi- 
spectral sensors. The current availability of multispectral imagery taken 
over the same geographical area on a time sequential basis (e.g. a month 
apart) is quite limited, however. The objective of this portion of the study 
was to identify and acquire time sequential multispectral imagery in as many 
spectral bands and covering as many terrain types as possible. 

The NASA Earth Resources Research Data Facility was used extensively for this 
imagery search. Multispectral imagery taken from the NASA Earth Observations 
Aircraft Program were found to be the most applicable. The frames of multi- 
spectral imagery selected for digitization are identified in Table 1-1 
(references 1, 2, and 3). This imagery was selected because of good mission 
overlap in three spectral bands . 
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TABLE 1-1 

SELECTED FRAMES OF MULTISPECTRAL IMAGERY 


Mission 

No. 


Site 

No. 


General 

Location 


Roll 

No. 


Frame 

No. 


Filter 

No. 


Spectral 

Band 


Terrain 

Category 


131 


210 


Imperial Valley, 


14 


26 


25A 


Red 



131 


210 


California 



16 


26 


89 B 


Green 


Farmland 


131 


210 



- 


15 


26 


58 


IR 




131 


29 


City of 



14 


54 


25A 


Red 




131 


29 


Phoenix 



16 


54 


89B 


Green 


Urban 


131 


29 





15 


54 


58 


IR 




131 


29 


Vicinity of 



14 


97 


25A 


Red 




131 


29 


Phoenix, Ariz 


16 


97 


89 B 


Green 


Mountain 


131 


29 





15 


97 


58 


IR 




131 


217 


Vicinity of 



4 


45 


25A 


Red 




131 


217 


Atlanta, Georgia 


5 


45 


89 B 


Green 


Forest 


131 


217 





6 


45 


58 


IR 




129 


210 


Imperial Valley, 


7 


16 


25A 


Red 




129 


210 


California 



8 


16 


89 B 


Green 


Farmland 


129 


210 





9 


16 


58 


IR 




129 


29 


City of 



7 


51 


25A 


Red 




129 


29 


Phoenix 



8 


51 


89 B 


Green 


Urban 


129 


29 





9 


51 


58 


IR 




129 


29 -1 

Vicinity of 



7 


57 


25A 


Red 




129 


29 


Phoenix, Arizona 


8 


57 


89B 


Green 


Mountain 


129 


29 





9 


57 


58 


IR 




158 


217 


Vicinity of 



31 


48 


25 


Red 




158 


217 


Atlanta, Georgia 


34 


48 


89B 


Green 


Forest 


158 


217 





32 


48 


58 


IR 
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Mission 129 was flown on 21 May 1970 at all three sites. Mission 131 was 
flown on 16 June 1970 for sites 210 and 29, and 8 June 1970 for site 217. 
Mission 158 was flown on 5 March 1970. Additional missions are available at 
monthly time increments for sites 29 and 210 which were not included in this 
study. 

The terrain types selected from the available imagery are Farmland, Urban, 
Mountain, and Forest. Of these, the farmland terrain imagery was the easiest 
to identify. The forest terrain imagery was the most difficult. The selected 
frame for the forest terrain could, in fact, be labeled "miscellaneous”, as 
much open land was included as well as wooded areas. 

Imagery representative of the total visible spectrum (color) and the blue 
spectrum were not available for these missions. 

1.3 RESOLUTION REQUIREMENTS 

Image ground resolutions of 10 meters and 50 meters were selected as being 
typical requirements for advanced earth observation platforms. The imagery 
acquired from the NASA Earth Resources Research Data Facility was 70 mm. 
black and white negative transparencies. Typically, this imagery was photo- 
graphed from an altitude of approximately 18,288 meters (60,000 ft.). This 
section discusses how the 70 mm. transparencies were photo-processed and 
digitized to meet the approximate 10 meter and 50 meter resolution require- 
ments . 

Without specific ground truth for the selected frames of imagery, it is im- 
possible to determine image resolution with any degree of accuracy. By 
using a ground resolution formula (reference 4) however, which relates camera 
lens resolution, film resolution, and scale, a ground resolution upper bound 
for a photographic image can be determined (depending on the film type used) 
as follows: 
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(f) (R^) (3.28084) vx x ' 

where : 

V = ground resolution in meters (per line pair) . 
f = focal length (mm.) 

= lens/film resolution (line pair/mm.). 

Using this formula, the upper resolution bound on the acquired 70 mm. imagery 
can then be determined and is noted in Table 1-2. 


TABLE 1-2 

RESOLUTION UPPER BOUND AS A FUNCTION 
OF FILM TYPE 


x 


Film Type 


r lf 

Line Pair/mm 


Ground 

Resolution 

(meters) 


I-. 1 2424 1 

80 


5.72 


2. 2448 


80 


5.72 


3. 2402 


100 


4.57 


4. 3400 


160 


2.86 
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It is apparent (from Table 1-2) that the acquired multispectral imagery has 
greater resolution than is required for this study. 

The MCAIR Digital Image Processing System (DIPS) (See Figure 1-2) can be used 
to filter the multispectral imagery as it scans and digitizes the image data. 
For equal vertical and horizontal resolution, -a square non-overlapping scanning 
lattice such as is illustrated in Figure 1-1, would be utilized. 
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FIGURE 1-1 

EFFECTIVE DIGITAL IMAGE PROCESSING SYSTEM 
SCANNING LATTICE FOR 8x8 ARRAY 
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In this scanning mode, the DIPS has a vertical and horizontal resolution 
of 13.1 lp/mm. If the DIPS were used to scan and digitze the 70 mm. imagery 
without any photo processing, the resulting ground resolution upper bound 
can be determined using equation (1-1) and setting = 13.1 lp/mm. Thus: 

V DIPS = 34.87 meters (1-2) 

This result implies that, if the DIPS were airborne at 18,288 meters (60,000 
feet) , scanning a fixed field of view, the maximum achievable resolution of the 
digital image data would be 34.87 meters. To increase or decrease the achievable 
ground resolution requires that the altitude parameter in Equation (1-1) be 
modified. This modification is most easily accomplished by suitable photo 
processing of the imagery before scanning. 

1.3.1 Photo Processing 

Thus, the filtering characteristic of the DIPS in the scanning mode, together 
with photo processing of the imagery before scanning, can be effectively 
used to provide digitized image data at the required resolutions. 
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a) To achieve 10 meter (32.8 feet) ground resolution: Photographi- 

cally enlarge the 70 mm imagery to 34.87/10 = 3.5 times the 
original size before scanning. 

b) To achieve 50 meter (164 feet) ground resolution: Photographically 

reduce the 70 mm imagery to 34.87/50 = .7 times the original 

size before scanning. 

Extreme care was used in the photoprocessing of the selected multispectral 
image frames. In addition to the enlargement or reduction of the imagery 
to achieve resolution requirements, scale differences (due to minor airplane 
altitude variations) were minimized in the photo processing step. 

1*3.2 Image Scanning and Digitization 

1. 3.2.1 Digital Image Processing System 

The MCAIR Digital Image Processing System (DIPS) was used to scan, sample, 
quantize and store the selected frames of multispectral imagery on computer 
compatible magnetic tape. For this study, the image intensity was quantized 
to 6-bit accuracy (0 - 63 levels). Image intensity was scanned instead of 
the more commonly used image density because of the image statistics require- 
ments for the image sampling study . (Section 3). 


The DIPS is a high resolution object plane image scanner and recording system. 
A block diagram of this system is illustrated in Figure 1-2. It has the 
capability of converting a photographic transparency into a rectangular or 
square sampling lattice of discrete image samples as defined by: 

a = ° ° a = °° 

I f s (x l’ X 2 ) = Kl ^ 6(x l " aa l ) ' l x 2 1 6(x 2 ” aa 2 ) * f < x i» x 2 ) 

a = — 00 a = — 00 

(1-3) 

Where: fCx^x^ is the continuous image function 
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Image Scanning Mode) 



FIGURE 1-2 

BLOCK DIAGRAM OF DIGITAL IMAGE PROCESSING SYSTEM 

GP72-0450-5 

The sampled function coin P ose ^ °f an arra Y of delta (6) functions 

spaced at intervals of width in the x^ direction and a 2 in the direction 
as shown in Figure 1-3. A value proportional to the intensity of the image 
area under each 6 function is recorded on magnetic tape by the DIPS. Any real 
world image scanner does not sample an image with ideal 5 functions. The 
measured sampling aperture of the DIPS is illustrated in two dimensional form 
in Figure 1-4. As can be seen, the sampling aperture is approximately 
Gaussian and the aperture width at half amplitude is 21 angstroms. 

For the image statistics to be representative of the image data, a» linear in- 
tensity response is required of the scanner. The intensity response of the 
DIPS in the scanning mode was measured and the required corrections are listed 
in Table 1-3. This intensity compensation was accomplished by table look up 
in the computer before the image statistics were determined * (see Figures 
1-7 and 1-25). As can be seen in Table 1-3, the intensity compensation causes 
a reduction in the number of quantized levels available. This could have been 
avoided by changing to 7-bit (128 levels) or 8-bit (256 levels) quantization 
for the image statistics computer processing. This was not done for this 
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FIGURE 1-3 

IDEAL SAMPLED IMAGE FUNCTION 


Legend: 

Dips Scanning Aperture 
O Gaussian Scanning Aperture f (x)] 


. JM i , x 2 

f(x) = r~ exp 

V 27r 2 0.737 






0.80 MS I 
61 angstromsj 
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FIGURE 1-4 

PLOT OF DIPS SCANNING APERTURE 
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TABLE 1-3 

DIGITAL IMAGE PROCESSINGJSYSTEM INTENSITY COMPENSATION 


DIPS | 

Corrected 
Quantization 
Level for 
Linear Intensity 
Variation 



DIPS! 

Corrected 
Quantization 
Level for 
Linear Intensity 
Variation 

Scanned 

Quantization 

Level 


Density 

Value 


Scanned 

Quantization 

Level 


Density 

Value 


0 


■9 


0 



32 


KISH 


29 


1 


mm 


2 



33 




30 


2 


1.05 


4 



34 




31 
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0.94 
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35 




32 
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0.84 
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36 


0.17 


33 
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0.78 
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37 




33 
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0.72 
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38 
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0.67 


11 



39 
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0.108 


40 


14 
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study because of time limitations and computer . programming problems. The 
intensity compensation in Table 1-3 was used for the study. 

1. 3.2.2 Generation of Multispectral Image Data 

One set of image data with an upper bound resolution of 50 meters was digitized 
from a portion of each frame of multispectral imagery listed in Table 1-1. 

This digital image library includes image data from three missions, and is 
representative of four terrain types in three spectral bands. In addition, 
one set of image data with an upper bound resolution of 10 meters was digitized 
for the red spectrum, all terrain categories, and all missions. DIPS reproduced 

versions of these images are displayed in Figures 1-8 through 1-23 and Figures 
1~26 through 1-41. It should be noted that the reproduced image quality shown 

has been degraded by both the DIPS in transparency reconstruction and photo 
processing for publication purposes. 

The DIPS was used exclusively for all image scanning and digitization. Digital 
image size, for all data images, was selected to be 128 x 128 samples. This 
choice was a compromise between the desirability of using a large image size so 
the image statistics are representative, and keeping the image size small to 
provide economical Fast Fourier Transform (FFT) computer processing. 

Rotational error between sequential flight imagery was minimized by mounting 
a 50X microscope on the DIPS to aid in image prescanning registration. This 
rotational error is estimated to be within +4 angstroms which is << one 
scanning aperture width. 

1.4 GENERAL MULTISPECTRAL IMAGE STATISTICS 

Figure 1-5 is a block diagram of the procedure used for generating the general 
multispectral image statistics. It is not clear what the best statistical 
descriptors for multispectral imagery are. Those selected for this study 
include the image intensity mean, variance, amplitude histogram, power spectral 
density function and autocorrelation function. It was anticipated that the 
image mean, variance, and autocorrelation statistics would be used in Section 
2 and the power spectral density used in Section 3. 
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FIGURE 1-5 

BLOCK DIAGRAM FOR GENERATING 
GENERAL IMAGE STATISTICS 


1.4.1 Development of Computer Program 

The power spectral density <j)(co) of an image can be defined as follows: 

+ («) - X(u>) # X*(w) (1-4) 

where X(a>) is the two dimensional Fourier transform of the 

spatial image function. 

X*(w) is the complex conjugate of the Fourier transform 
The autocorrelation function R^, x^) of an image can then be defined: 

R(x 15 x 2 ) = F _1 | X(u) X*(u) | (1-5) 

where : F ^ denotes the inverse two dimensional Fourier transform. 

Direct application of Equations (1-4) and (1-5) to numerical computation using 
the discrete Fourier transform must be performed with care (references 5, 6, 7) 
Cyclical correlation and "leakage" computational problems associated with the 
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discrete Fourier transform must be minimized before computational efficiencies 
using the Fast Fourier Transform (FFT) can be realized. According to Anuta 
(Reference 5) , the cyclical correlation problem is eliminated by padding the 
digitized image array with zeros (as shown in Figure 1-6) before transform 
computer processing. This padding greatly increases the size of the image 
data array and, hence, the computational cost. The leakage problem was 
minimized by first subtracting the image mean intensity from the data, and 
then multiplying the outer 10% of the 128 x 128 image array by the cosine bell 
function (References 6 and 7) . 





FIGURE 1-6 

’ — 1 

r \ 

ILLUSTRATION OF ZERO PADDING USING FFT 



j 
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j/ The computer program for generating the general image statistics is outlined 
' in Figure 1-7. The power spectral density function and autocorrelation 
. function arrays were determined using Equations (1-4) and (1-5) with zero 
padding of the image data and "leakage” minimization being accomplished. It 
J should be noted that $ (&) and R(x^,x^) are now 256 x 256 arrays because of 
^ the zero padding operation. 

1.4.2 Image Statistics 

General image statistics were determined for the 50 meter resolution and 10 
meter resolution digital images of Mission 131. These . statistics .are illus- 
trated in Figures 1-8 through 1-23. The three dimensional illustrations for 
the power spectral density and autocorrelation function arrays are plots of 
every eighth value in the array. 


T 
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FIGURE 1-7 

GENERAL IMAGE STATISTICS COMPUTER FLOW DIAGRAM 
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1 28 x 128 I mage Array 


32° 08’30"N. Latitude 
115° 28'00" W. Longitude 

Image Mean: 24.016 

Image Variance: 127,66 



Centered At: 


a. Reference Image 



b. Image Intensity Histogram 


FIGURE 1-8 

IMAGE STATISTICS (FARMLAND TERRAIN - RED SPECTRUM - 
MISSION 131 - 50 METER RESOLUTION) 
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d. Autocorrelation Function 


FIGURE 1-8 (Continued)) 

IMAGE STATISTICS (FARMLAND TERRAIN - RED SPECTRUM 
MISSION 131 - 50 METER RESOLUTION) 
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128 x 128 Image Array 


a. Reference Image 



Quantization Level 
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b. Image Intensity Histogram 


FIGURE 1-9 

IMAGE STATISTICS {FARMLAND TERRAIN - GREEN SPECTRUM - 
MISSION 131 - 50 METER RESOLUTION} 


MCOOW/VELL jtlRCJTAFT 

1-16 







Per Cent 


REPORT MDC A1740 



128 x 128 Image Array 
a. Reference Image 



b. Image Intensity Histogram 


FIGURE MO 

IMAGE STATISTICS (FARMLAND TERRAIN ■ IR SPECTRUM - 
MISSION 131 - 50 METER RESOLUTION) 
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Centered At: 33° 28'45" N. Latitude 

112° 04'45" W. Longitude 


< 


Image Mean: 17J965 

image Variance: 22.281 1 
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128 x 128 Image Array 

a. Reference Image 
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FIGURE 1-11 

IMAGE STATISTICS {URBAN TERRAIN - RED SPECTRUM - 
MISSION 131 - 50 METER RESOLUTION) 


Quantization Level 
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b. Image Intensity Histogram 


MCDOWWEU. AIRCRAFT 


1-20 


03 


d. Autocorrelation Function 


FIGURE 1-11 (Continued)[ 

IMAGE STATISTICS (URBAN TERRAIN - RED SPECTRUM 
MISSION 131 - 50 METER RESOLUTION) 
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128 x 128 Image Array 
a. Reference Image 
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b. Image Intensity Histogram 


FIGURE 1-12 

IMAGE STATISTICS (URBAN TERRAIN - GREEN SPECTRUM - 
MISSION 131 - 50 METER RESOLUTION) 
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128 x 128 Image Array 


Centered At: 


Image Mean: 


33° 28'45" N. Latitude 
112° 04’45" W. Longitude 

8.58197 


Image Variance: 12 r 639l 


a. Reference Image 



FIGURE 1-13 

IMAGE STATISTICS (URBAN TERRAIN - IR SPECTRUM - 
MISSION 131 - 50 METER RESOLUTION) 
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128 x 128 Image Array 
a. Reference Image 



b. Image Intensity Histogram 
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FIGURE 1-14 

IMAGE STATISTICS (MOUNTAIN TERRAIN - RED SPECTRUM - 
MISSION 131 - 50 METER RESOLUTION) 
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c. Power Spectra! Density Function 
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d. Autocorrelation Function 
FIGURE 1-14 (Continued)] 

IMAGE STATISTICS (MOUNTAIN TERRAIN - RED SPECTRUM 
MISSION 131 - 50 METER RESOLUTION) 
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128 x 128 Image Array 
a. Reference Image 


Centered At: 

Image Mean: 
Image Variance: 


33° 30'45" N. Latitude 
1 1 1° 31 '45" W. Longitude 

11.801 

5. 1 3853 



b. Image Intensity Histogram 


FIGURE 1-15 

IMAGE STATISTICS (MOUNTAIN TERRAIN - GREEN SPECTRUM - 
MISSION 131 - 50 METER RESOLUTION) 
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d. Autocorrelation Function 
FIGURE 1-15 (Continued)f 

IMAGE STATISTICS (MOUNTAIN TERRAIN - GREEN SPECTRUM - 
MISSION 131 - 50 METER RESOLUTION) 
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Centered At: 

Image Mean: 
Image Variance: 


33° 30'45" N. Latitude 
111° 31'45" W. Longitude 

13.55 

11.0946 


1 28 x 128 Image Array 


a. Reference Image 
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b. Image Intensity Histogram 


FIGURE 1-16 

IMAGE STATISTICS (MOUNTAIN TERRAIN - IR SPECTRUM - 
MISSION 131 - 50 METER RESOLUTION) 
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FIGURE 1-16 (Continued)| 

IMAGE STATISTICS (MOUNTAIN TERRAIN - IR SPECTRUM - 
MISSION 131 - 50 METER RESOLUTION)! 
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1 28 x 128 Image Array 
a. Reference Image 



GP72 0450-19 

b. Image Intensity Histogram 


FIGURE 1-17 

IMAGE STATISTICS (FOREST TERRAIN - RED SPECTRUM - 
MISSION 131 - 50 METER RESOLUTION) 
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d. Autocorrelation Function 
FIGURE 1-17 ( Continued)j 

IMAGE STATISTICS (FOREST TERRAIN - RED SPECTRUM - 
MISSION 131 - 50 METER RESOLUTION) 
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128 x 128 Image Array 


a* Reference Image 



Quantization Level 


GP72 0450-18 


b. Image Intensity Histogram 


FIGURE 1-18 

IMAGE STATISTICS (FOREST TERRAIN - GREEN SPECTRUM - 
MISSION 131 - 50 METER RESOLUTION) 
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128 x 128 Image Array 


a. Reference Image 



b. Image Intensity Histogram 


FIGURE 1-19 

IMAGE STATISTICS {FOREST TERRAIN - IR SPECTRUM - 
MISSION 131 - 50 METER RESOLUTION) 
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d. Autocorrelation Function 


FIGURE 1-19 (Continued)^ 

IMAGE STATISTICS (FOREST TERRAIN - IR SPECTRUM - 
MISSION 131 - 50 METER RESOLUTION) 
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128 x 128 I mage Array 


a. Reference Image 



b. Image Intensity Histogram 


FIGURE 1-20 

IMAGE STATISTICS (FARMLAND TERRAIN - RED SPECTRUM - 
MISSION 131 ■ 10 METER RESOLUTION) 
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128 x 128 Image Array 


Image Mean: 12.7518 

Image Variance: 8.16483 


a. Reference Image 
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FIGURE 1-21 b 


FIGURE 1-21 

IMAGE STATISTICS (URBAN TERRAIN - RED SPECTRUM - 
MISSION 131 - 10 METER RESOLUTION} 
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d. Autocorrelation Function 


FIGURE 1-21 (Continued) 

IMAGE STATISTICS (URBAN TERRAIN - RED SPECTRUM - 
MISSION 131 - IOMETER RESOLUTION) 
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128 x 128 Image Array 
a. Reference Image 



Quantization Level 
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b. Image Intensity Histogram 


FIGURE 1-22 

IMAGE STATISTICS (MOUNTAIN TERRAIN - RED SPECTRUM - 
MISSION 131 - 10 METER RESOLUTION) 



MCOOWJVfLt AIRCRAFT 

1 - 1*2 





REPORT MDC A1740 



Image Mean: 64 

Image Variance: 10.5 


128 x 128 Image Array 


a. Reference Image 



b. Image Intensity Histogram 


FIGURE 1-23 

IMAGE STATISTICS (FOREST TERRAIN * RED SPECTRUM - 

MISSION 131 - IOMETER RESOLUTION) 
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d. Autocorrelation Function 


FIGURE 1-23 (Continued)[ 

IMAGE STATISTICS (FOREST TERRAIN - RED SPECTRUM - 
MISSION 131 - 10 METER RESOLUTION) 
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1.5 DIFFERENCE IMAGE STATISTICS 


When properly generated, a difference image contains any change information 
between two images of the same geographical area at two separate observation 
times. This change information is a result of the following: 

a) Different sun angles (shadows) 

b) Different viewing angles 

c) Seasonal changes 

d) Sampling 

e) Random sensing and data processing system phenomenon 

Figure 1-24 is a block diagram of the procedure used for generating the 
difference image statistics. 



FIGURE 1-24 

BLOCK DIAGRAM FOR GENERATING DIFFERENCE IMAGE STATISTICS 
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1.5.1 Development of Computer Program 


The flow diagram for the computer program used to generate the difference 
image and its relevant statistics is illustrated in Figure 1-25. 
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The program input consists of two digital image arrays from magnetic tape. 

The (x^, x^) displacement between the two images is determined by utilizing 
the cross-correlation coefficient image registration algorithm. Computation- 
ally, this was accomplished using two dimension FFT techniques. Figures 
l-26b through l-41b l are three dimensional plots of the center (32 x 32) 
points for this image registration estimator. (In some cases, the image 
alignment point is not shown in these plots.) Following determination of the 
(x , x ). displacement, the data image is forced to register with the refer- 
ence image. The overlapped portions of the two images are then identified 
and the mean intensity value of the overlapped images are made equal. The 
difference image is then formed and the relevant difference image statistics 
printed out. 

1.5.2 Image Statistics 

Difference image statistics were determined for the 50 meter resolution and 
10 meter resolution digital images scanned in Section 1.3. 2. 2. The differ- 
ence images are defined and identified in Table 1-4. The relevant statistics 
determined were the difference image mean, variance, and the quantized prob- 
ability distribution discussed more fully in Section 2.3.1. 2. These 
statistics are illustrated in Figures l-26c and d through Figures l-41c 
and d. 
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TABLE 1-4 

DIFFERENCE IMAGE IDENTIFICATION 


DIFFERENCE IMAGE 
(Reference Image - Data Image) 


Figure No 

• 

Figure No. 

l-8a 

Minus 

l-26a 

l-9a 

ii 

l-27a 

l-10a 

ii 

l-28a 

1-lla 

ii 

l-29a 

l-12a 

ti 

l-30a 

l-13a 

it 

1-3 la 

l-14a 

ii 

l-32a 

l-15a 

ii 

l-33a 

l-16a 

ii 

l-34a 

1-1 7a 

it 

l-35a 

1-1 8a 

ii 

l-36a 

l-19a 

it 

l-37a 

l-20a 

n 

l-38a 

l-21a 

ii 

1-3 9a 

l-22a 

it 

l-40a 

l-23a 

it 

l-41a 
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128 x 128 Image Array 


a. Data Image 


c. Difference Image Statistics 
(Reference Image ■ Data Image) 



b. Image Registration Estimator 


FIGURE 1-26 

DIFFERENCE IMAGE STATISTICS (FARMLAND TERRAIN - RED SPECTRUM - 
MISSION 129 - 50 METER RESOLUTION) 
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Quantization Level - x 10 


d. Difference Image Histogram 


FIGURE T26 (Continued)! 

DIFFERENCE IMAGE STATISTICS (FARMLAND TERRAIN - RED SPECTRUM - 
MISSION 129 - 50 METER RESOLUTION) 
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128 x 128 Image Array 
a. Data Image 


Difference Image Mean: 0,326 

Difference Image Variance: 17.1 


c. Difference Image Statistics 
(Reference Image - Data Image} 



b. Image Registration Estimator 


FIGURE 1-27 

DIFFERENCE IMAGE STATISTICS (FARMLAND TERRAIN - GREEN SPECTRUM - 
MISSION 129 • 50 METER RESOLUTION} 
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128 x 128 Image Array 
a. Data Image 


Difference image Mean: —0,095 

Difference image Variance: 40.3 


c. Difference Image Statistics 
(Reference Image ■ Data Image) 




b. Image Registration Estimator 


FIGURE 1-28 

DIFFERENCE IMAGE STATISTICS (FARMLAND TERRAIN - IR SPECTRUM - 

MISSION 129 - 50 METER RESOLUTION) 
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4.0 

3.9 * 



d. Difference Image Histogram 


> FIGURE 1-28 (Continued)/ 

DIFFERENCE IMAGE STATISTICS (FARMLAND TERRAIN - IR SPECTRUM - 
MISSION 129 - 50 METER RESOLUTION) 
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128 x 128 Image Array 
a. Data Image 


Difference Image Mean: 0,087 

Difference Image Variance: 6.07 


c. Difference Image Statistics 
(Reference Image - Data Image) 



b. Image Registration Estimator 


FIGURE 1-29 

DIFFERENCE IMAGE STATISTICS (URBAN TERRAIN - RED SPECTRUM - 
MISSION 129 - 50 METER RESOLUTION 
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Quantization Level - x 10 


d. Difference Image Histogram 


FIGURE 1-29 (Continued)! 

DIFFERENCE IMAGE STATISTICS (URBAN TERRAIN - RED SPECTRUM - 
MISSION 129 - 50 METER RESOLUTION 
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128 x 128 Image Array 


Difference Image Mean: 0,48 

Difference Image Variance: 17,2 



a. Data Image 


c. Difference Image Statistics 
(Reference Image - Data Image) 



b. Image Registration Estimator 


FIGURE 1-30 

DIFFERENCE IMAGE STATISTICS (URBAN TERRAIN - GREEN SPECTRUM 
MISSION 129 - 50 METER RESOLUTION) 
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Quantization Level - x 10 


d. Difference Image Histogram 


FIGURE 1-30 (Continued)! 

DIFFERENCE IMAGE STATISTICS (URBAN TERRAIN - GREEN SPECTRUM - 
MISSION 129 - 50 METER RESOLUTION) 
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128 x 128 Image Array 
a. Data Image 


Difference Image Mean: —0.12 

Difference Image Variance: 16.3 


c. Difference Image Statistics 
(Reference Image - Data Image} 



b. Image Registration Estimator 


FIGURE 1-31 

DIFFERENCE IMAGE STATISTICS (URBAN TERRAIN - IR SPECTRUM - 
MISSION 129 - 50 METER RESOLUTION} 
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d. Difference Image Histogram 


FIGURE 1-31 (Continued)) \ 

DIFFERENCE IMAGE STATISTICS (URBAN TERRAIN - IR SPECTRUM 
MISSION 129 - 50 METER RESOLUTION) 
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128 x T 28 Image Array 


a. Data Image 


Difference Image Mean: —0.21 

Difference Image Variance: 4.15 


c. Difference Image Statistics 
(Reference Image* Data Image) 




b. Image Registration Estimator 


FIGURE 1-32 

DIFFERENCE IMAGE STATISTICS (MOUNTAIN TERRAIN - RED SPECTRUM - 
MISSION 129 - 50 METER RESOLUTION) 
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d. Difference Image Histogram 


FIGURE 1-32 (Continued)[ 

DIFFERENCE IMAGE STATISTICS (MOUNTAIN TERRAIN - RED SPECTRUM - 
MISSION 129 - 50 METER RESOLUTION) 
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128 x 128 Image Array 
a. Data Image 


Difference Image Mean: 0.459 

Difference Image Variance: 7.13 


c. Difference Image Statistics 
(Reference Image - Data Image} 



FIGURE 1-33 

DIFFERENCE IMAGE STATISTICS (MOUNTAIN TERRAIN - GREEN SPECTRUM - 
MISSION 129 - 50 METER RESOLUTION) 
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d. Difference Image Histogram 


FIGURE 1-33 (Continued)[ 

DIFFERENCE IMAGE STATISTICS (MOUNTAIN TERRAIN - GREEN SPECTRUM . 

MISSION 129 - 50 METER RESOLUTION) I 
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128 x 128 Image Array 


Difference Image Mean: 0.37 

Difference Image Variance: 6.03 


a* Data Image 


c, Difference Image Statistics 
(Reference Image* Data Image) 



b. Image Registration Estimator 


FIGURE 1 34 
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d. Difference Image Histogram 


FIGURE 1-34 (Continued)[ 
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128 x 128 Image Array 
a. Data Image 


Difference Image Mean: 0.44 

Difference Image Variance: 52,2 


c. Difference Image Statistics 
(Reference Image- Data Image) 



b. image Registration Estimator 


FIGURE 1-35 

DIFFERENCE IMAGE STATISTICS (FOREST TERRAIN - RED SPECTRUM - 
MISSION 158 - 50 METER RESOLUTION) 
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d. Difference Image Histogram " ’ 

FIGURE 1-35 (Continued)[ 

DIFFERENCE IMAGE STATISTICS (FOREST TERRAIN - RED SPECTRUM 
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1 28 x 128 1 mage Array 


a. Data Image 


c. Difference Image Statistics 
{Reference Image* Data Image) 



b. Image Registration Estimator 


FIGURE 1-36 

DIFFERENCE IMAGE STATISTICS {FOREST TERRAIN - GREEN SPECTRUM - 
MISSION 158 - 50 METER RESOLUTION) 
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Quantization Level - x 10 


d. Difference Image Histogram — ' 

FIGURE 1-36 (Continued)| 

DIFFERENCE IMAGE STATISTICS (FOREST TERRAIN - GREEN SPECTRUM - 
MISSION 158 - 50 METER RESOLUTION) 
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Difference Image Mean: 0.082 

Difference Image Variance: 30.5 


128 x 128 Image Array 


a. Data Image 


c. Difference Image Statistics 
(Reference Image - Data Image) 



b. Image Registration Estimator 


FIGURE 1-37 

DIFFERENCE IMAGE STATISTICS (FOREST TERRAIN - IR SPECTRUM - 
MISSION 158 - 50 METER RESOLUTION) 
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d. Difference Image Histogram 
FIGURE 1-37 ( Continued)! 

DIFFERENCE IMAGE STATISTICS (FOREST TERRAIN - IR SPECTRUM - 
MISSION 158 - 50 METER RESOLUTION) 
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128 x 128 Image Array 


a. Data Image 


Difference Image Mean: 0.211 

Difference Image Variance: 57.6 


c. Difference Image Statistics 
(Reference Image- Data Image) 



b. Image Registration Estimator 


FIGURE 1-38 

DIFFERENCE IMAGE STATISTICS (FARMLAND TERRAIN - RED SPECTRUM - 
MISSION 129 - 10 METER RESOLUTION) 
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d. Difference Image Histogram 
FIGURE 1-38 (Continued)| 

DIFFERENCE IMAGE STATISTICS (FARMLAND TERRAIN - RED SPECTRUM - 
MISSION 129 - 10 METER RESOLUTION) 
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128x 1281 mage Array 



Difference Image Mean: 0.359 

Difference Image Variance: 3,59 


a. Data Image 


c. Difference Image Statistics 
(Reference Image - Data Image) 



b. Image Registration Estimator 


FIGURE 1-39 

DIFFERENCE IMAGE STATISTICS (URBAN TERRAIN - RED SPECTRUM - 
MISSION 129 - 10 METER RESOLUTION) 
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d. Difference Image Histogram 


FIGURE 1-39 (Continued)) 

DIFFERENCE IMAGE STATISTICS (URBAN TERRAIN - RED SPECTRUM - 
MISSION 129 - 10 METER RESOLUTION) 
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128 x 128 Image Array 
a. Data Image 


Difference Image Mean: 0.072 

Difference Image Variance: 8.16 


c. Difference Image Statistics 
(Reference Image - Data Image) 



b. Image Registration Estimator 


FIGURE 1-40 

DIFFERENCE IMAGE STATISTICS (MOUNTAIN TERRAIN - RED SPECTRUM - 
MISSION 129 - 10 METER RESOLUTION) 
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Quantization Level - x 10 


d. Difference Image Histogram 


FIGURE 1-40 (Continued)! 
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128 x 128 Image Array 
a. Data Image 



Difference Image Mean: -0,5 

Difference Image Variance: 8 


c. Difference Image Statistics 
(Reference Image ■ Data Image) 



b. Image Registration Estimator 


FIGURE 1-41 

DIFFERENCE IMAGE STATISTICS (FOREST TERRAIN - RED SPECTRUM - 
MISSION 158 - 10 METER RESOLUTION) 
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SECTION 2 

IMAGE CORRELATION STUDY 


2.1 INTRODUCTION 

In this study, digital techniques for the two-dimensional spatial registration 
of multispectral imagery are developed and evaluated. * As the imagery from 
successive flights shown in Section 1 has demonstrated, there are two basic 
reasons why spatial registration is a serious problem for potential users of 
multispectral imagery. 

1) The Fields of View (FOV) are not precisely aligned. 

2) The images are sensed at different times from different sensor 
locations. 

These factors cause translation, rotation, scale, geometric distortion, and 
brightness/contrast differences between images to be registered. User appli- 
cations of earth resources imagery require accurate spatial alignment so 
that multidimensional analysis can be performed on image data from different 
spectral bands, or from successive overflights, or both. 

This study considers the translation problem associated with digital image 
registration and develops a means for comparing commonly used digitial 
algorithms. Using suitably defined constraints, an optimum and five suboptimum 
position location (registration) techniques are developed and evaluated. Both 
Gaussian and image derived statistics are used to compare the various tech- 
niques as functions of radial error (in image sample distances) and computa- 
tional cost. 
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2.2 IMAGE REGISTRATION GEOMETRY 


This section describes the image geometry, nomenclature and assumptions used 
to examine the image translational registration problem. 

Two digital picture element arrays are defined as follows: 


S = Reference Image (m x M) 


(2-1) 


d = Data Image (n x N) 


(2-2) 


where d is considered to be a subset of S and (n x N) is smaller than (m x M) . 
These arrays are illustrated in Figure 2.1 where the (n x N) data image is 
referenced to its equivalent position on the reference image, that is, the 
potential alignment point (£,L) is the position on the reference image of the 
(1, 1) position on the data image. 


t i 



Mf 



(M)| 

Tj 

n| 

L 



1 

- — n| — -1 


a. Reference Image 

(s e+ i - 1, l + j -i) 


b. Data Image 

(d u> 


FIGURE 2-1 ! 

IMAGE REGISTRATION GEOMETRY 


"\ 

J 
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In the usual image registration search procedure, the data image is systemati 
cally stepped in two dimensions around the reference image. This procedure 
implies that the position location alignment point (£° , L°) must lie within 
the range' 1 i _< M— N+l and 1 <_ j _< m-n+1 on the reference image 

The difference image, formed by subtracting S from d, is sometimes called 
the "noise" or "change detected" image* These image changes are caused by 
varying sun angles, different viewing angles, seasonal changes and other 
random phenomena associated with the sensing system. For this study, the 
difference image £ is defined as follows: 


e, ,a,L) 

x J J 



S £+i-l, L+j-P 


(2-3) 


where represents additive sample-to-sample independent noise. This 

assumption makes it possible to derive the optemum position location 
technique disscussed in Section 2*3.1. 


2.3 POSITION LOCATION ANALYSIS 

In the following sections, an optimum position location algorithm based on 
minimum mean square radial error (SRE) estimation is derived and its com- 
putational complexity discussed. In addition, five suboptimum position 
location algorithms commonly used for image registration are defined together 
with their computational requirements . 


2.3.1 The Optimum Position Location Technique 

The conditional probability of alignment on the reference image S (assuming no 
angular misalignment) shown in Figure 2-1 is given by Equation (2-4). 

P(£,L) = CP u,L) n n P {c> a, L) } (2-4) 

i=l j=l 4. i,j 
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Where : C = normalization parameter. 

a priori probability of alignment for any given point in S. 


C. . (A,L)= defined in Equation (2-3). 
1 9 J 


n N 

FI II P { £ . . (^,L) }= joint probability of the difference image £ • 
i=l j-1 4 


For this study, no a priori knowledge about the alignment point (£°,L°) is 
known, therefore, P (£,L) is assumed 1/ (m-n+1) (M=N+1) . The probability 
of alignment P(£,L) can be determined for each allowable point within S 
for any given difference image £. Allowable points (£,L) consist of all 
points within the set f i 1 ^ 1 m-n+1, 1 <iX _< M-N+1 } . Assume that the 
alignment. point (&° , L°) is within the allowable point set. The normalization 
parameter C is then specified as: ' 


n 

o 


m-n+1 M-N+1 

X=1 ■ '■ L=1 


n N 

n n 

i-i j-i 




-1 


(2-5) 


Since: 


m-n+1 


z 

■ 0 = 1 


M-N+1 

y: p(a,d a i 

L=1 


(2-6) 


2. 3. 1.1 Mean Square Radial Error 


The mean square radial error (MS RE) estimate for the alignment point is 
determined for that position of (£,L) which minimizes: 


r (Jl 1 ,L 1 ) = 


m-n+1 M-N+1 

E E 

“ L=1 


P(£,L) [(i^SL) 2 + (L 1 ~L) 2 J 


(2-7) 


Where: r (£,L) = mean square radial error. 

P(£,L) = defined in Equation (2-4). 

(A^,L^) = alignment point being examined. 
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The implementation of this optimum position location estimator for image 

registration applications requires a tractable expression for P {?. . (£,L)} 

b 1 9 J 

to solve for P(£,L) in Equation (2-4). In the following section, an approach 
for implementation of this estimator is described. 


2 . 3 . 1 . 2 Implementation 


One attractive approach for implementation of the optimum position 

location estimator is to determine P {?. . (£,L)} experimentally using 

b 1 >3 

histograms of real difference image data and then perform a minimization of 


r2(£,L) (Equation 2-7) using the digital computer. Unfortunately, when this 
approach is implemented using reference images of useful size, the digital 
storage and computational cost of using this estimator becomes prohibitive. 

An alternate approach and the one determining P^ { } in a straightforward 
manner used in this study is to form the following equivalent expression. Let 


n 

all i 


n 

all j 


p r (c, .a,u) = n 

■ all q 


b (£,L) 

q 


(2-8) 


Where: b (£,L) = the number of difference image samples that lie in the 

qth quantization level. 


P^ = the probability that a difference image sample is in the 
qth quantization level. 

Using the right half of Equation (2-8) instead of the left to determine 
P(£,L) allows a substantial savings in both image statistic storage and 
estimator computations. A typical quantized difference image probability 
P function is illustrated in Figure 2-2. 

q 

For any given difference image then, (assuming the constraints of Sections 
2.2 and 2.3.1) the optimum position location algorithm can be digitally 
implemented for image registration using the following steps: 

a) Form P histogram for given £ . . (£,L). 

q hj 

b) Solve for T1 , 11 P {?. . (&,L)} using Equation (2-8). 

alii all j C 1,J 
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c) Solve for normalization constant C using Equation (2-5). 

d) Determine P‘(£,L) using Equation (2-4). 


e) Solve for the alignment point (£*, L*) where r (£,L) is minimum 
:hat the only a ] 
the difference image. 


Note that the only a priori image statistic required is the P^ histogram of 





2. 3.1. 3 Estimator 

The optimum position location estimator was described mathematically in 
Equation (2-7). For image registration, the alignment point is that value 
of (£,L) where r2(£,L) is a minimum. Figure 2-3 is a 3-D plot illustrating 
the behavior of this estimator for a specific difference image. For the 
example of Figure 2-3, n=N=5, m=M=50, and a difference image with a Gaussian 
Pq histogram (a = 10) was assumed. 
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2 . 3 . 1 . 4 Computations 

The optimum position location algorithm requires a specific number of compu- 
tations for any given reference and data image size. These computations have 
been listed parametrically in Table 2-1. In this tabulation, a subtraction 
and division were made equivalent to an addition and multiplication respec- 
tively. 
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TABLE 2-1 

OPTIMUM POSITION LOCATION ALGORITHM COMPUTATIONS 


Parameter 


Additions 


Multiplications 


n N / 

jj, P f fu (S - L 


nN, V (£, L) 


(nN-1), V (£, L) 

2. C, Store Value 


(m-n+1) (M-N + 1) 


0 

3. P (fi, L). Store Each Value 


nN, V (£, L) 


nN, V (£, L) , 

4. r 2 (£, L) 


4(m-n + 1) (M-N + 1) -1, 
V (£, L) 


3(m— n + 1) (M— N + 1), 
V (£, L) 

5. All r 2 (£, L) 


4(m— n + 1) 2 (M-N + 1) 2 + 
+ nNj(m-n + 1)(M-N + 1) 

1 

1 

3(m— n + 1) 2 (M-N + 1) 2 
+ nN (m— n + 1) (M— N + 1) 
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-U- 

2.3,2 v « 1 Correlation Techniques 

The v-law correlation algorithms are based on the minimization of the 
absolute difference between S and d. For example, if S and d are identical 
(at alignment), the sum of the absolute values of the difference image 
samples are zero. The difference image at other potential alignment points 
would have non- zero samples yielding the separating property of the v-law 
correlation algorithm. Alignment for the general v-law correlation algorithm 
occurs for that point (£*,L*) where: 


n 

N 



i=l 

l d i,j ” S SL+±-l, 

r 

L+j -1 

(2-9) 

is a minimum. 

For v=l, this becomes: 



n 

N 



«i- z 

S l d i,j ” S £+i-l, 

L+j-1 1 

(2-10) 

i=l 

j-1 
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2. 3. 2.1 Estimator 

The v-1 position location estimator described in Equation (2-10) can best be 
illustrated by a specific examples. Figure 2-4 is a 3-D plot of the estimator 
using a Gaussian difference image (a/A=10), n=N=5, m=M=50 (See Section 2.4). 



2 . 3 . 2 . 2 Computations 


For any given set of S and d images, a finite number of computations are 
required using the v=l position location estimator. No multiplications are 
required, only additions. These are: 

Number of Additions 
For one (£,L) 

Total Number of Additions 
For All Allowable (£,L) 


| = (2nN-l) (m-n+1) (M-N+l) (2-12) 


= 2nN-l 


( 2 - 11 ) 
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2*3.2 v = 2 Correlation Technique 

Another attractive suboptiriium position location technique involves using 
Equation. (2-9) with v=2. 

n N 

C 2 = S ^ d i» j ” S i+i-l, L+j-ll (2_ 

i=l j=l 

2. 3. 3.1 Estimator 

Alignment occurs for that point (£*,L*) where C£ is a minimum. A 3-D plot 
of a representative v=2 position location estimator is shown in Figure 2-5. 
As with the previous estimators, this example assumes a Gaussian difference 
image (o/A=10), n=N=5 , m=M=50 (See Section 2.4). 



EXAMPLE OFi 


Li , 

FIGURE 2-5 1 

v=2 POSITION LOCATION ESTIMATOR 
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2. 3. 3. 2 Computations' 

For any given set of S and d images, a finite number of computations 
are required using the v=2 position location estimator. These are: 

Total Number of Additions ) 

For All Allowable (£,L) j = (2nN-1) (m-n+1) (M-N+l) (2-14) 

Total Number of Multiplications j 

For All Allowable (£,L) j (m-n+1) (M-N+l) (2-15) 




2.3.4 Maximum Likelihood Estimation Technique 

The maximum likelihood estimator is, perhaps, the most utilized decision 
algorithm. Applying this estimation technique to the image registration 
problem, the maximum likelihood position location algorithm would compute 
the probability of alignment at all allowable points in S. The potential 
alignment point (,£*,L*) having the maximum probability value would be 
selected as the alignment point. The expression for computing these align- 
ment probabilities was given previously in Section 2.3.1, Equation (2-4). 

n N. 

pa,L) = c n n p r u. ,a,L)> \2-uj 

i=l j=i 5 1,J 


When the difference image 4 is Gaussian with zero mean, the maximum 
likelihood estimator becomes identical to the v=2 law position location 
algorithm. This can be shown as follows: 

Let: 


n 

n 

i=l 


N 

n 

j=i 


P H t i,j (t,L)} 



(2-17) 


Where: C = normalization constant discussed in Equation (2-5), 


c 2 - 


n 


N 


n2 


^ d i»j ~ S £+i-l, 

Li=l j=l 


L+j-l) 


(2-18) 
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Using maximum likelihood estimation, the point (£*,L*) would be selected as 

the alignment point which maximizes P(£,L) in Equation (2-16). This procedure 

2 

is equivalent to minimizing ? as defined in Equation (2-18) which is seen 
to correspond with Equation (2-13) for the v=2 algorithm. Thus, for Gaussian 
zero mean £, the maximum likelihood and v=2 are equivalent position location 
algorithms. This equivalence will, be noted in Section 2.4. 

2. 3.4.1 Estimator 

Due to the equivalence of the maximum likelihood and v=2 algorithms when £ 
is a zero mean Gaussian difference image, Figure 2-5 illustrates a represen- 
tative maximun likelihood position location estimator. 

2. 3. 4. 2 Computations 

The finite number of computations required by the maximum likelihood position 
location technique is determined from Table 2-1 to be; 

Total Number of J 

Additions For all = (nN+1) (m-n+1) (M-N+l) (2-19) 

Allowable (£,L) ' 

Total Number of \ 

Multiplications for = (nN+1) (m-n+1) (M-N+l) (2-20) 

All Allowable (£,L) ) 

r' . 

2.3.5 Cross-Covariance Function Techniques 

The cross covariance position location estimator is defined as: 
n N 

C xy " ^2 ^2 (d i,j -d) (S £+i-l, L+j-l"^,^ (2-21) 

i=l j=l 

Where: d = mean Value of d. . 

_ 1,3 

S £,L = Mean Value of S £+i-l, L+j-1 ’ summed over all i and j 


’ f 
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The algorithm alignment point estimate occurs for that point where 

C xy has its maximum value. 

2. 3.5.1 Estimator 

A 3 — D plot of a representative cross- covariance position location estimator is 
shown in Figure 2-7. This example assumed a Gaussian difference image 
(ct/A= 10) , n=N=5, m=M=50 (See Section 2.4) 



FIGURE_2-6| | 

EXAMPLE OF CROSS-COVARIANCE POSITION LOCATION ESTIMATOR i 


2. 3* 5. 2 Computations 


The expression for the cross-variance position location estimator can be re- 
written as follows: 


xy 


^ d i,j S £+i-l, L+j-l) 

V=1 j-1 / 


nN d S £,L (2 ~ 22) 
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Equation (2-22) can be computed with fewer additions and multiplications 
than Equation (2-21) and will be used to determine computational require- 
ments. (A subtraction is equivalent to an addition.) 


Total Number of Additions 
For All Allowable (£,L) 


(m-n+1) (M-N+l) (2nN-l) + nN-1 


(2-23) 


Total Number of \ 

Multiplications for >= (m-n+1) (M-N+l) (nN+4) + 1 
All Allowable (&,L) 


2.3.6 Cross Correlation Coefficient Technique 

The most widely used digital position location algorithm is that of the 
cross-correlation coefficient (Reference 5 ) . This estimator is defined m 

Equation (2-25) in terms of two dimensional image samples. 



Where: C is defined by Equation (2-22) 

xy 

d = Mean value of d 

1 J 3 

= Mean Value of S £+i-l, L+j-1 

Alignment occurs at that point (A*,L*) where Pjcy has its maximum value. 


2»3.6»1 Estimator 

A 3-D plot of a representative cross-correlation coefficient position 
location estimator is shown in Figure. 2-8. This example assumed a 
Gaussian difference image (a/ A = 10), n=N=5, and m=M-50 (See 
Section 2.4) . 
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FIGURE 2-7 

EXAMPLE OF CROSS-CORRELATION COEFFICIENT POSITION 
LOCATION ESTIMATOR 


2. 3.6.2 Computations 

Using conventional numerical methods, the computational requirements for the 
cross-correlation coefficient position location technique are listed as 
follows. A subtraction and division are equivalent to an addition and multi- 
plication, respectively. 

Total Number of I 

Additions For j = (4nN-l) (m-N+1) (M-N+l)+2nN (2-26) 

All Allowable Points J 

Total Number of I 

Multiplications for j = (2nN+6) (m-n+1) (M-N+1 )+17N+3 (2-27) 

All Allowable Points I 
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In addition, Fast Fourier Transform (FFT) computational methods can be 
utilized for correlation. Such usage however, must be exercised with care. 
Aliasing errors require the d image be padded with zeros to form an m x M 
image. An additional constraint is that m and M must be powers of 2. For 
determining computational requirements using the FFT, it is convenient to 
assume n=N and m=M. Then: 

Total Number of I T w _ o 

> = 6M Z Log 2 M + 2M Z (2-28) 

Additions j 


Total Number of j 

Complex Multiplications J * 3 « 2 L ^2 « + + »< 2 +l < 2 - 29 > 

2.3.7 . Position Location Optimization 

This section discusses various parameters which affect both the accuracy and 
efficiency of position location techniques. 


2 , 3 . 7 . 1 Preprocessing 

Preprocessing is defined as .any refinement of the image data S and d which 
will reduce the system error in the £ image. This includes any photo- 
processing and image normalization which make £ representative only of terrain 
’’change detected” or difference information. 


Preprocessing is very important to the optimum, v=l, v=2, and maximum like- 
lihood position location techniques. These estimators, basically, all search 
for minimum values of the £ image. At alignment, any error due either to 
dynamic range limitation or intensity mean value offset will seriously 
degrade performance of these position location techniques. The cross- 
correlation position location techniques have some normalization inherent 
in the algorithms which reduce their susceptibility to these errors. 

2 . 3 . 7 . 2 Quantization 

The optimization of image quantization for image registration applications 
has a twofold purpose: 1) reduce the error probability and 2) reduce the 

computational complexity. Examination of each position location algorithm 
discussed in the preceding sections will show that each is based on a set of 
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sums of random variables. The sum accuracy for each algorithm is a direct 
function of the quantization range for that random variable. For example, 
the estimator for the optimum position location technique (Equation 2-7) is 
seen to be a set of sums whose accuracy depends on the quantization error 
present in the ? image (Equation 2-3). 


The error probability for a sum of random variables can be approximated by 
invoking the central limit htheorem: this error probability is: 

1 n n’ 7 



(2-30) 


where: d, . 

1,3 


= unquantized data image 


d. *. 
(d. *. 


V«> 


= quantized data image 

y v 

- d. .) = quantization error e. . (d) is 

1 » J x,j 

uniformly distributed 

= probability that the sum is outside the confidence interval 

[-a, a] . 


Consider the normalized random variable: 


y = 



(2-31) 


Which has a limiting distribution that is normal with zero mean and unity 
variance. In other words, as "nN" becomes larger, the probability density 
function for y is approximately a zero mean unity variance normal. Thus, 
the probability that |y| is greater than some 6 is: 
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( n 

N / 

) 


H|y|>«) = p< 1 Y. 


m 

(2-32) 

{ i=i 

2 

3=1 




Wh ere: A_ = variance of uniformly distributed quantization error 

12 (reference 8 ) 

2 

A nN = Variance of a sum of uniformly distributed independent 
12 quantization error random variables. 

A- = quantization interval 


P(|y|>6)- 





(2-33) 


Equation (2-33) and (2-31) are now equivalent and are plotted in Figure 2-8. 

It is obvious from Figure 2-8 that, for a fixed a, increasingly fine quanti- 
zation will decrease the error probability. This implies that increasingly 
fine quantization should provide image registration with increasing position 
location accuracy. A trade off of registration accuracy as a function of 
quantization probably exists for each algorithm, however, a digital simulation 
would be needed to provide quantitative evaluation. 

2.3. 7.3 Data Image Size and Shape 


In general, an optimum data image size and shape should exist for each of 
the position location techniques. Intuitively, those position location 
estimations that depend on the correlation present between two images (such 
as the and algorithms) should benefit from using maximum sample 

compactness in the sampling lattice. This would imply that a rhombic 
sampling lattice (as discussed in Section 3.4.3) should be used for these 
algorithms. It is felt that the improvement over using a square sampling 
lattice would only be slight, however. Quantitive data is lacking concerning 
the trade off of data image size and shape as a function of position location 
error performance. In the digital simulations of Section 2-4 and 2-5, an 
attempt is made to determine if square or circular data images should be 
used. 
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FIGURE 2-8| 

ERROR PROBABILITY AS A FUNCTION OF QUANTIZATION INTERVAL] 
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2. 3. 7. 4 Computations 

In this section, a comparison of computational requirements for the position 
location techniques is presented. This computational comparison assumes 
that "brute force" numerical techniques are used as specified in the previous 

sections. Table 2-2 is a collection of the required computation formulae. 

TABLE 2-2 

COMPUTATION FORMULAE 


Position 

Location 

Algorithm 


Computation Operations] 


Equivalent Additions] 

Equivalent Multiplications] 

1. Optimum 



4A.j 2 A 2 2 + nN A 1 A 2 


3 Ay 2 A 2 2 + nN A 1 A 2 


2. v= 1| 



(2 nN -1) A 1 A 2 




3. v= 2) 



(2nN — 1 ) A 1 A 2 


n(\l A2 


4. Maximum Likelihood 


(nN + 1) A 1 A 2 


(nN + 1) A 1 A 2 


I 5 - C xy| . " 



(2nN — 1 ) A 1 A 2 + nN — 1 


(nN + 4) A, A 2 + 1 


6. p xy (Conventional^ 


(4nN — 1) A 1 A 2 + 2nN 


(2nN + 6) A 1 A 2 + nN + 3 

1 

7 - /V FFTI l 



6M 2 Log 2 M + 2M 2 


3M 2 Log 2 M + A 2 2 + M 2 + 

l| 

where: = 

A2 = 

m— n + 1 
M-N + 1 
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For a straightforward comparison, the computation "cost" for a multiplication 
has been set at 3 times that required for an addition. Other assumptions 
include the use of square reference and data images. The curves of Figure 
2- 9 ' demonstrate this comparison as a function of data image size for a 
reference image m=M=50. 

It is apparent from these curves that the optimum position location tech- 
nique requires the largest number of computations and the v=l technique 
requires the least number. v = 2 and C^are almost identical. 

To compare p using the FFT with the other algorithms, it is necessary to 
use a reference i m age size of a power of 2. Table 2-3 is such a comparison 
using m=M=64 and varying n=N=8, 16 and 32. Because of the zero padding 
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FIGURE 2-9 

COMPUTATION COMPARISON 

(Reference Image m = M - 50)i 
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requirement for the data image, the number of computations for using 
the FFT is basically due to the size of the reference image and does not 
vary appreciably with data image size. For large data images (e.g. n=N=32) , 
the p estimator using the FFT is seem to require far less computations 
than the other position location techniques. 


TABLE 2-3 

COMPUTATION COMPARISON 

(Reference Image m = M = 64) 


n = N| 

Number of 
Samples 
in Data 
Image 

Number of Additional Operations Required) 

Optimum) 

^=1| 

n 

Maximum 

Likelihood 

m 

■1 

n 

8 

16 

32 

1 

64 

256 

1024 


0.138 x 10 9 
0.774 x 10 8 
0.199 x 10 8 

0.413 x 10 6 
0.123 x 10 7 
0.223 x 10 7 


0.845 x 10 6 
0.247 x 10 7 
0.446 x 10 7 

0.213 x 10 7 
0.619 x 10 7 
0.101 x 10 8 


0.108 x 10 7 
0.31 x 10 7 
0.559 x 10 7 
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2. 3. 7. 5 Thresholds 

Various thresholds can be employed in certain position location algorithms 
to reduce computational effort. The adaptive threshold shown in Figure 2-10 
could be employed for those position location estimators that look for a 
minimum value (such as the optimum, v=l, and v =2 techniques). An additional 
requirement is that the estimators should be well behaved functions (e.g., 
compare the estimators of Figures 2-3, 2-4 and 2-5 with those of Figures 2-6 
and 2-7). The adaptive threshold is based on the comparison of each new 
estimator value against the previously computed minimum value. Computations 
for a position location estimator would terminate when the threshold is 
exceeded. For the case where the estimator has a final value which is less 
than the threshold, the adaptive threshold would assume the minimum computed 
value. The computational savings achieved by use of the adaptive threshold 
can only be determined through simulation. The simulations in Section 2-4 
and 2-5 used such a threshold and results are reported there. 
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In addition, variable threshold schemes can be used to reduce computational 
effort. Reference [9_ is an example of this kind of approach. It should be 
noted however, that the use of a variable threshold has the potential for 
increasing the probability of error. 



2. 3* 7* 6 Search Scanning Philosphy 

The optimum starting point for any of the position location estimators is 
the center of the allowable points in S. There are several searching al- 
gorithms that can be utilized to speed the searching process. For example 

a coarse/fine algorithm can be easily emplemented for the C and p esti- 

xy xy 

mators. Such searching algorithms were not investigated in this study since 
their use would lead to an increase in the probability of position location 
error. 

2.4 POSITION LOCATION USING GAUSSIAN DIFFERENCE IMAGES 

As real difference image statistics did not become available until the latter 
part of this study, a digital simulation was developed to compare the position 
location techniques using Gaussian computer generated difference images. Trade 
offs of data image size and shape were examined as a function of radial error 
performance and estimator value for difference image noise levels of 'a/ A = 1, 
2.5> 5, and 10. Reference image size was selected to be 50 x 50 and is the 
center 50 x 50 subimage of the reference image illustrated in Figure l-8a. 
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Both square and circular data images were evaluated, 
the lattice used for the d. . images when n=N=15 and 


Figure 2-11 illustrates 
A (diameter) =15 . 


2.4.1 Development of Computer Program 

The flow diagram for the position location simulation is presented in Figure 
2-12. This computer program determines position location performance as a 
function of radial error (in image sample distances) for 64 level and 2 level 
(hardlimited) quantized images. 


The Gaussian P distribution was machine generated using the NOISE SUBROUTINE v 
^ ^ , 
to form difference images having noise levels of a/ A = 1, 2.5, 5, and 10 

quantization intervals. Statistically, these ? images can be compared to 
real world difference images having variances o 2 = 1, 6.25, 25, and 100. 

The data image d is formed by adding the difference image £ to the reference 
image S (shown in Figure 1.8a). Subjectively, d resembles the images illus- 
trated in Figure 2-13 for a / A = 1, 2.5, 5, and 10. To reduce computation 
requirements, only the center 50 x 50 image samples were used in this simu- 
lation. 

It should be noted that, for statistically reliable results using this 
simulation, a large number of different noise images should be run at each 
specified a / A level. Such large numbers of simulation runs, however, were be- 
yond the scope of the current effort. For the results presented in the 
following three sections, only one such machine generated noise image was 
utilized. These results demonstrate that the computer program is operational 
and can be used to perform trade off studies for the indicated image registra- 
tion parameters. In this simulation, the maximim likelihood estimator was 
not included due to its similarity with the V = 2 estimator as discussed in 
Section 2.3.4. 

L _ - 

~ i ...... tC z -- ^ 

2.4.2 Simulation Results 


2. 4. 2.1 Radial Error 


Position Location Radial Error is defined as: 

Radial Error = ^(£°-£*) 2 + (L°-L*) 2 

(£°, LP) is the alignment point 
(£*, L*) is the alignment point selected by the estimator 

f 
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a. Square Lattice (n = N = 15) 

Total Samples = 225 

o 

o o o o o o o 
ooooooooo 
ooooooooooo 
ooooooooooooo 
ooooooooooooo 
ooooooooooooo 
ooooooooooooooo 
ooooooooooooo 
ooooooooooooo 
ooooooooooooo 
ooooooooooo 
ooooooooo 

O O O O O O O ! 

o 

b. Circular Lattice (A = 15) 

Total Samples = 149 

FIGURE 2-11 I 

COMPARISON OF SQUARE AND CIRCULAR DATA IMAGE SHAPE G P72-0450-42 1 
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FIGURE 2-12 \ 

FLOW DIAGRAM OF POSITION LOCATION DIGITAL SIMULATION 
USING GAUSSIAN DIFFERENCE IMAGES 
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FIGURE 2-12 \(Continued)| 

FLOW DIAGRAM OF POSITION LOCATION DIGITAL SIMULATION 
USING GAUSSIAN DIFFERENCE IMAGES 
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c. a ! A = 5 



d. at A= 10 



FIGURE 2-13 

TYPICAL NOISY DATA IMAGES (n = N = 128) 
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The computer program outlined in Figure 2-12 determines the radial error for 
the optimum and suboptimum position location algorithms for a specified 
difference image. Using a single computer generated difference image, data 
image parameters of size and shape were examined as functions of a/ A, where: 



a_ - 1 sigma value of g (2-35) 

A quantization internal 


For this simulation, A = 1; therefore, a and a/ A were varied from 1, 2.5, 
5 , and 10 . 


Radial error simulation results are illustrated in Tables 2-4, 2-5, 2-6, and 
2-7. Although it is difficult to draw conclusions from using only a single 
difference image, a few general observations can be made. 

a) The optimum position location estimator did not always provide 
the least radial error. 

b) Except for the a / A = 1 case, radial error did not seem to be a 
function of a/ A. 

c) The C position location algorithm is the suboptimum technique 
most likely to provide registration etror. 

d) In many cases, 2-level quantization provided similar radial 
error performance to 64-level -quantization (but not all 
cases! e.g., see of A = 10). 

e) A data image size of less than 150 samples appears to be 
adequate for image registration. 

These comments apply to the Gaussian difference image simulation only. 

2; 4:2; 2 Estimator Values 

The estimator values at alignment are given in Tables 2-8, 2-9, 2-10, and 
2-11. This information is of value in establishing typical thresholds for deter- 
mining image registration versus no registration. For those cases in the tables 
where 0 radial error resulted, the following observations can be made: 

a) An estimator value of less than .5 appears to provide image 
registration with 0 radial error for the optimum algorithm. 

Probably an estimator value of less than .1 would be a better 
choice. 
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TABLE 2-4 


RADIAL ERROR COMPARI 


SQUARE DATA IMAGE 


POSITION 

LOCATION 

ALGORITHM 

n = N 

Data 

Image 

Size 

RADIAL ERROR 
(In Sample Distances) 

64 Level 

2 Level 

r 2 (£,L) 
v = 1 
v = 2 
C 

p Xy 
. xy 

5 

25 

25 

25 

25 . 

25 

25 

0 

0 

0 

0 

11.7 

0 

0 

0 

0 

0 

22 

7 

r 2 (Jt,L) 

15 

225 

0 ‘ 

0 

• v = 1 


225 

0 

0 

v - 2 


225 

0 

0 

C 


225 

0 

0 

xy 

P 


225 

0 

13 

xy 


225 

0 

13 

r 2 (*,L) 

25 

625 

0 

0 

v — 1 


625 

0 

0 

v = 2 


625 

0 

0 

C 


625 

0 

0 

xy 


625 

0 

0 

P xy 


625 

0 

0 


CIRCULAR DATA IMAGE 


Diameter 


Data 

Image 

Size 


RADIAL ERROR 
(In Sample Distances) 


64 Level 


2 Level 


17.26 

22.09 

22.09 

22.09 

0 

7 
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POSITION 

LOCATION 

ALGORITHM 


TABLE 

RADIAL ERROR COMP 


SQUARE DATA IMAGE 


r (*,L) 
v = 1 
v = 2 
C 

xy 


rV.L) 

v “ 1 

v = 2 

C 

xy 

P 

xy 


r(*,L) 
v = 1 
v = 2 
C 

xy 

P xy 


5 

25 

19.24 

18.44 


25 

19.42 

19.24 


25 

19.24 

19.24 


25 ’ 

19.24 

19.42 


25 

15.81 

19.42 

■ 

25 ‘ 

19.24 

7 

15 

225 

0 

0 


225 

0 

0 


225 

0 

0 


225 

0 

0 


225 

0 

13 


225 

0 

13 

25 

625 

0 

0 


625 

00 

0 


625 

0 

0 


625 

0 

0 


625 

0 

0 


625 

0 

o- 
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TABLE 2-6 

RADIAL ERROR COMPARISON FOR a/ A = 5 




SQUARE 

DATA IMAGE 


CIRCULAR 

DATA IMAGE 

POSITION 


Data 

RADIAL 

ERROR 


t a 

RADIAL 

ERROR 

LOCATION 

ALGORITHM 

n = N 

Image 

(In Sample Distances) I 

A 

Image 

(In Sample 

Distances) 



Size 




Size 








Diameter 






64 Level 

2 Level 



64 Level 

2 Level 

r 2 (£,L) 

5 

25 

15.30 

5.39 

5 

13 

11.66 

4.47 

V) = 1 


25 

14.14 

0 


13 

22.36 

25.46 

v = 2 


25 

21.1 

0 


13 

18.68 

25.46 

c 


25 

21.1 

16.03 


13 

18.68 

27. 66 

n Xy 


25 

11.7 

14.87 


13 

22.47 

0 

P 


25 

7.28 

7 


13 

20.22 

9 

xy 









r 2 (l,L) 

15 

225 

0 

0 

15 

n 

0 

0 

v = 1 


225 

0 

0 


BB 

0 

0 

v = 2 


225 

0 

0 


BUB i 

0 

0 

C 


225 

0 

0 


149 

0 

0 

xy 

P 


225 

17.03 

13.15 


149 

10.05 

0 

xy 


225 

0 

13 


149 

0 

5 

r 2 (^,L) 

25 

625 

0 

0 

25 

441 

0 

0 

v = 1 


625 

0 

0 


441 

0 

0 

v = 2 


625 

0 

0 


441 

0 

0 

C 


625 

0 

0 


441 

0 

0 

xy 


625 

0 

0 


441 

0 

0 

P xy 


625 

0 

0 


441 

0 

0 



I 
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POSITION I 
LOCATION I n = N 
ALGORITHM I 


r (Z,L) 
v = 1 
v = 2 
C 

xy 


r (A,L) 
v = 1 
v = 2 


r (*,L) 
v = 1 
v = 2 
C 


Data 

Image 

Size 


DATA IMAGE 

RADIAL 
(In Sample 

ERROR 

Distances) 

64 Level 

2 Level 


CIRCULAR DATA IMAGE 


Data 

Image 

Size 


RADIAL ERROR 
(In Sample Distances) 


64 Level I 2 Level 


5 

25 

12 

11 


25 

12.4 

29. 7 


25 

13.3 

29.7 


25 

13.3 

15.26 


25 

23.02 

28.6 


25 

17.12 

7 

15 

225 

0 

1.41 


225 

0 

-4.12 . 


225 

0 

4.12 


225 

0 

1.41 


225 

17 

13v04 


225 

0 

0 

25 

625 

0 

1 


625 

0 

0 


625 

0 

0 


625 

0 

1 


625 

0 

0 


625 

0 

0 


13.34 

24.04 

13.34 

13.34 

23.77 

18.11 


6 . 

24. 
24. 
24.04 
0 
7 




12.17 

12.17 

12.17 

12.17 

0 

2 
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TABLE 2-8 

COMPARISON OF ESTIMATOR VALUES FOR ct/A = 1 




SQUARE 

DATA IMAGE 


CIRCULAR 

DATA IMAGE | 

POSITION 


Data 




Data 



LOCATION 

n = N 

Image 

ESTIMATOR VALUES 


Image 

ESTIMATOR VALUES 1 

ALGORITHM 


Size 




Si 7, p 






64 Level 

2 Level 



64 Level 

2 Level 

r 2 (£,L) 

5 

25 

.22xl0 _14 

.03 

5 

13 

.61 

144 

V = 1 


25 

29 

1 


13 

18 

0 

v = 2 


25 

61 

1 


13 

44 

0 

Q 


25 

1 

.999 


13 

.939 

.717 

xy 


25 

571 

6.2 


13 

244 

6 

P 

xy 


25 

.78 

.14 


13 

.792 

.926 

r 2 (»,L) 

15 

225 

0 

. 73xl0~ 16 

15 

mm 

0 

, -12 
.53x10 

v = i 


225 

221 

11 



142 

7 

v = 2 


225 

371 

11 



228 

7 

C 


225 

_ 1 

1 


mm- 

1 

1 

xy 

o 


225 

31000 

46 


149 

15600 

99 

xy 


225 

.994 

.627 


149 

.992 

.147 

r 2 (*,?0 

25 

625 

0 

0 


441. 

0 

.11x10 

V - 1 


625 

630 

'.17 

25 

441 

439 

13 - 

v = 2 


625 

1130 

17 


441 

777 

13. 

C 


625 

1 

1 


441 

1 

1 

xy 


625 

97200 

145 


441 

55300 

203 

P xy 


626 

.994 

.976 


441 

.993 

.116 
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TABLE 2-9 

COMPARISON OF ESTIMATOR VALUES FOR a/A = 2.5 




SQUARE 

DATA IMAGE 


CIRCULAR 

DATA IMAGE 

1 

POSITION 


Data 




Data 



LOCATION 

n = N 

Image 

j ESTIMATOR VALUES 


Image 

ESTIMATOR VALUES 1 

ALGORITHM 


Size 



Size 






64 Level 

2 Level 


64 Level 

2 Level 

r 2 (£,L) 

5 

25 

1.54 

18 

5 

13 

17.2 

275 

V = 1 


25 

58 

3 


13 

32 

1 

V = 2 


25 

234 

3 


13 

116 

1 

C 


25 

.685 

.835 


13 

.91 

.014 

xy 


25 

686 

4.8 


13 

314 

3 

P 

xy 


25 

.707 

.109 


13 

.819 

.79 

r 2 ( £,L) 

V = 1 

15 

225 

225 

0 

551 

.002 

27 

15 

149 
• 149 

0 

348 

.0028 

19 

X 

v = 2 


225 

2180 

27 


149 

1300 

19 



225 

1 

.999 


149 

1 

.999 

L« 

xy 


225 

31100 

43 


149 

15500 

91 

P 


225 

.967 

.569 


149 

.958 

.139 

xy 








r 2 ( W 

25 

625 

0 

.195xl0~ 38 

25 

441 

0 

. 31x10“ 

v = 1 

625 

1620 

49 

441 

1130 

38 

v = 2 


625 

6860 

49 


441 

4710 

38 

C 


625 

1 

1 


441 

1 

1 

xy 


625 

97400 

129 


441 

54600 

189 

P xy 


' 625 

.967 

.547 


441 

.959 

.114 
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TABLE 2-10 

COMPARISON OF ESTIMATOR VALUES FOR a! A = 5 




SQUARE 

DATA IMAGE 


CIRCULAR 

DATA IMAGE | 

POSITION 


Data 




Data 



LOCATION 

n a N 

Image 

| ESTIMATOR VALUES 

A 

Image 

ESTIMATOR VALUES 1 

ALGORITHM 


Size 




Si z p 






64 Level 

2 Level 

Diameter 


64 Level 

2 Level 

rV,L) 

5 

25 

98 

169 

5 

13 

165 

351 

v = l 


25 

124 

6 


13 

52 

2 

v = 2 


25 

971 

6 


13 

319 

2 

c 


25 

.216 

.22 


13 

.0964 

.0121 

xy 

A 


25 

1800 

4.68 


13 

671 

3 

K 

xy 


25 

.713 

.14 


13 

.8128 

.1080 

r 2 (£,L) 

15 

225 

-22 

.36x10 

. 8x10 ^ 

15 

149 

.3647x10 ^ 

5.123 

v = 1 


225 

1160 

33 


149 

766 

30 

v = 2 


225 

9100 

33 


149 

5942 

30 

C 


225 

1 

.999 


149 

1 

.9609 

xy 

A 


225 

31100 

38.3 


149 

15505 

80 

K 

xy 


225 

.878 

.509 


149 

.84 

.11845 

rV ,L ) 

25 

625 

0 

.13x10 -2P 

25 

441 

-45 

.596x10 

. 23xl0 -17 

v = 1 


625 

3290 

87 


441 

2305 

69 

v = 2 


625 

26800 

87 


441 

18211 

69 

C 


625 

1 

1 


441 

1 

1 

xy 


625 

98100 

111 


441 

56563 

171 

P xy 



.888 

.44 


441 

.87128 

.4458 
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TABLE 2-11 

COMPARISON OF ESTIMATOR VALUES FOR a ! A = 10 



* 

3 




SQUARE 

DATA IMAGE 


CIRCULAR 

DATA IMAGE | 

POSITION 


Data 




Data 



LOCATION 

n = N 

Image 

ESTIMATOR VALUES 

A 

Image 

ESTIMATOR VALUES 1 

ALGORITHM 


Size 




Size 






64 Level 

2 Level 

Diameter 


64 Level 

2 Level 

r 2 («,,L) 

5 

25 

80.7 

208 

5 

13 

176.5 

313 

v = i 


25 

201 

8 


13 

103 . 

2 

v = 2 


25 

2470 

8 


13 

154 

2 

C 


25 

.075 

.00313 


13 

/ 02446 

.006158 

xy 

A 


25 

1310 

3.6 


13 

867 

3 

y 

xy 


25 

.547 

.11 


13 

' .696. 

.926 

r 2 (£,L) 

15 

225 

.0235 

1.519 

15 

mm 

.02855 

.09654 

v = 1 


225 

2130 

57 


m 

1395 

41 

v = 2 


225 

31500 

57 


■9 

19725 

41 

C 


225 

.994 

.88900 


149 

.9926 

.9786 

xy 

o 


225 

30100 

32.15 


149 

15097 

60 

K 

xy 


225 

.687 

.284 


149 

.6133 

.1020 

r 2 (^,L) 


625 

.73xl0 _16 

.153 


441 

181 

. 4.9xl0 -1 

v = 1 

25 

625 

6150 

165 

25 

441 

4281 

118 

v = 2 


625 

95600 

165 


441 

65981 

118 

C 


625 

1 

.85 


441 

1 

.99972 

xy 


625 

101,000 

74.3 


441 

56297 

142 

P xy 


625 

.726 

.317 


441 

.68 

.105 
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b) There does not seem to be a well-defined image registration 
reject threshold for the -v = 1 and v = 2 algorithms. 

c) An estimator value of greater than 31,000 appears to provide 
image registration with 0 radial error for the cross-covariance 
algorithm. 

d) A well-defined image registration reject threshold was not 
demonstrated for the cross-correlation coefficient algorithm. 

Further runs would probably indicate an estimator value in the 
.8 to .9 range. 

Several additional simulation runs are needed to develop statistically 
reliable results. 

2. 4. 2.3 Computational Savings 

Adaptive thresholds were implemented in the simulation for the optimum, v = 1, 
and v = 2 position location estimators. These thresholds reduce the number of 
computations required for alignment without increasing the probability of error. 
Using results from the simulation for the single difference image, averaged 
computational savings were determined based on the numerical techniques dis- 
cussed in Section 2. 3. 7. 4. On the average, the optimum position location 
algorithm required only 23% (Q = 64) and 26% (Q = 2) of the equivalent 
addition computations indicated in Figure 2-10 and Table 2-3. The v = 1 
algorithm required 31% (for Q= 64 and Q = 2), and the v = 2 algorithm required 
19% (Q = 64) and 31% (Q = 2). It should be noted (from Table 2-3) that use of 
the adaptive threshold makes the v, = 1 position location estimator competitive 
with the (FFT) estimator on a computational basis. 

2.5 POSITION LOCATION USING MULTI SPECTRAL IMAGE STATISTICS 

This section examines estimator position location performance using difference 
image statistics taken from real multispectral imagery. These statistics 
(illustrated in Section 1.5) became available for use late in the study. 

2.5.1 Development of Computer Program 

The position location digital simulation of Section 2.4 can be used to accept 
statistics from real difference images. The difference image histograms 
determined in Section 1.5.2 for the selected multispectral imagery can be 


MCDO/V/VfLI. >Wf?CMFr 



REPORT MDC A1740 


entered directly into the computer program of Figure 2-13. It should be empha- 
sized that, for statistically reliable results, several noise images should be 
generated for each P^ distribution. Limitations of scope and time, however, 
allowed only one such noise image to be generated for the simulation results 
reported in the next several sections. 

2.5.2 Simulation Results 

As in Section 2.4, trade offs of algorithm performance as a function of data 
image size and shape were accomplished. In addition, one run was made to deter- 
mine algorithm radial error as a function of terrain category, image spectrum, 
and image resolution. 

2. 5. 2.1 Radial Error As Function of Data Image Size and Shape 

Circular and square data images were evaluated using the P histogram illus- 
trated in Figure l-26d. The variance of this difference image (Figure l-26c) 
is 63.1 and is approximately comparable to the Gaussian a/ A = 10 noise 
simulation. The reference image was taken to be the center 50 x 50 samples of 
Figure l-26a. A comparison of radial error performance for the various 
position location algorithms is given in Tables 2-12 and 2-13. 

Again, it is difficult to draw any substantive conclusions from the use of 
only one noise image. The radial errors recorded for the estimators are not 
unreasonable if the difference image is considered to be "worst case". It 
can be seen from examination of the difference image variances (Figures l-26c 
through l-41c) , that this indeed is the situation. Comparison can be made 
between these results and the a/A => 10 Gaussian noise simulation of Section 

2.4.1 (Table 2-7). Again, the optimum position location estimator did not 
always provide the least amount of registration error. In general (but not 
always) , the radial error was greater for Q = 2 quantization than for Q = 64 
for all algorithms except C . For C , radial error became less than that 
for Q = 64 with increasing data image size. 

2. 5. 2. 2 Radial Error As Function of Image Spectrum 

One run of the simulation was made to compare radial error performance as a 
function of image spectrum. Farmland terrain was assumed together with a 
square Nl = n = 7 data image. These results are given in Table 2-14. 
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TABLE 2-12 

RADIAL ERROR COMPARISON FOR SQUARE DATA IMAGE 


RADIAL POSITION LOCATION ERROR 


DATA 

IMAGE 

SIZE 


Optimum 


ikelihood! 


uantization = 64 Levels 


2.24 

10.05 

0 


17 

9.06 

19.1 

11 

12 


14.04 

10 

8.06 

0 

0 


25 

14.87 

2.83 

0 

1.41 


18.44 

14.21 

23.71 

1 

1 


1 

17.49 

2.24 

1 

1 


0 

17.12 

0 

2 

0 

0 

0 

13.89 


uantization = 2 Levels 


1 

0 

9.9 

10.05 

0 

13.89 

7 


14.14 
17.49 
15.03 

13.15 


14.14 
17.49 
15.03 

13.15 


19.24 

0 

24.17 

11.18 
0 

22.63 

7 


12.73 

11.4 

11.4 

13.04 


2 — 
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TABLE 2-13 

RADIAL ERROR COMPARISON FOR CIRCULAR DATA IMAGE' 



I RADIAL POSITION 

LOCATION ERROR 


DATA 

IMAGE 

SIZE 

Optimum 

r ~ 

|- 77 - 

Max 

Likelihooc 

P xy 

C 

xy 

Quantize 

Ltion = 64 

Levels 



13 

4.12 

14.04 

14.04 

13.04 

11.05 

14.21 

29 

5.1 

3.61 

3.61 

10.44 

18.44 

15.56 

49 

1.41 

1 

1 

10.44 

16.55 

15.56 

81 

4. 

13.04 

13.04 

21.4 

10.44 

21.10 

149 

25.46 

2.0 

14. 

15.03 

0 

13.04 

253 

19.85 

0 

0 

19.85 

0 

14.76 


Quantiza 

tion = 2 L 

evels 




13 

5.1 

14.04 

14.04 

14.04 

11. 

14.04 

29 

8.49 

3.61 

3.61 

3. 

11. 

21.47 

49 

6.71 

15.56 

15.56 

17. 

17. 

0 

81 

8.94 

16.28 

16.28 

1.41 

4. 

13.45 

149 

9.85 

13.6 

13.6 

8.54 

10. 

12.21 

253 

19.85 

0 

0 

19.85 

16. 

0 
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TABLE 2-14 

RADIAL ERROR COMPARISON AS FUNCTION OF IMAGE SPECTRUM 



| RADIAL POSITION LOCATION ERROR | 

Cmage 

Spectrum 

| Optimum 

v = 1 

v = 2 

Likelihood 

P 

xy 

c 

JSL 1 

| Quantization = 64 Levels 1 

Red 

Green 

IR 

8.06 

7.28 

5 

0 

14.14 

24.19 

0 

3.61 

24.19 

10 

15.03 

24.19 

25 

29 

24.19 

14.87 

28.43 

14.21 


Quantizal 

tion = 2 Le 

vels 




Red 

Green f 
IR 

9.9 

1.41 

8.6 

14.21 

23.09 

20.81 

14.21 

23.09 

20.81 

24.17 

18 

8.49 

8 

21 

17 

14.14 

10.20 

20.59 


Pq, histograms from Figures l-26d, l-27d, and l-28d were used. It should be 
noted (from Figures l-26c, l-27c, and l-28c) that the variance of the differ- 
ence images are 63.1 for the Red spectrum, 17.1 for the Green spectrum, and 

40.3 for the IR spectrum. 

2. 5. 2. 3 Radial Error As Function of Terrain Type 

One run of the simulation was made to compare radial error performance as a 
function of image terrain type. The Red spectrum was selected for registra- 
tion together with a square n = N = 7 data image. These results are presented 
in Table 2-15. P q histograms from Figures l-26d, l-29d, l-32d, and l-35d were 
used. The variances for these difference images are 63.1 for the Farmland 
terrain, 6.07 for the Urban terrain, 4.15 for the Mountain terrain, and 52.2 
for the Forest terrain. 

2. 5. 2. 4 Radial Error as Function of Image Resolution 

One run of the simulation was made to compare radial error performance as a 
function of image resolution. Again Farmland terrain was assumed together 
with a square n = N = 7 data image. The digital simulation, up to this point, 
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TABLE 2-15 

RADIAL ERROR COMPARISON AS FUNCTION OF TERRAIN TYPE 



Radial Position Location Error | 

Terrain 

Type 


Optimum| 

1| 

» = 2| 

Maximum 

Likelihood 

n 

^xyl 

C xy| 



Quantization = 64 Levels | 

Farmland 


8.06 


0 


0 


10 


25 


14.87 


Urban 



26.25 


27.66 


27.66 


26.25 


0 


0 


Mountain 


26.25 


0 


0 


26.25 


0 


16.28 


Forest 



1 


17.03 


8.06 


8.06 


8.06 


17 








Quantization 

= 2 Levels 





Farmland 


9.9 


14.21 


14.21 


24.17 


8 


14.14 


Urban 



1 


29.70 


29.70 


31.11 


0 


0 


Mountain 


12.04 


18.97 


18.97 


18.97 


0 


20.52 


Forest 



5 


19.34 


19.42 


14.42 


1 


19.42 



had used only 50 meter resolution image data. The 10 meter image data of 

i 

] Figure 1-38 is compared with the 50 meter image data of Figure 1-26 in Table 
I 2-16. The variances for the difference images are 63.1 for the 50 meter 
resolution, and 57.6 for the 10 meter resolution image data. 


TABLE 2-16 

RADIAL ERROR COMPARISON AS FUNCTION OF IMAGE RESOLUTION 


Image 

Resolution 

(Meters) 

Radial Position Location Error) 

Optimum! 

1! 

v = 2 \ 

Maximum 

Likelihood 

^xy 

1 

C xy| 

Quantization = 64 Levels) 

50 

10 

50 

10 

8.09 

12.81 

i 

0 

23.6 

i 

H 

10 

11 

25 

11 

i 

■ 

1 

Quantization =2 Levels! 

9.9 

5 

1 




24.17 

9.22 

8 

13 



1 
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2.6 CONCLUSIONS 

This study has derived a suitably defined optimum (on the basis of minimum 
squared radial error estimation) image registration algorithm. It has 
developed and provided the computational means for the quantitative compari- 
son of the optimum and five suboptimum image registration techniques using 
both machine generated and real multispectral difference image statistics. 
Computational requirements for the six image registration algorithms were 
determined and illustrated. Radial error performance of the algorithms was 
determined and illustrated using one run of machine generated and one run 
of real multispectral difference image statistics. Parametric trade offs 
of image quantization, and data image size and shape were demonstrated. 

On the basis of two statistical runs however, it is impossible to make 
major conclusions. Some general comments that can be made are: 

a) The optimum image registration algorithm appears to be a valid 

means for comparing the performance of suboptimum techniques. Further 
refinements can be made to improve its performance. 

b) On a computational basis, the v=l (using adaptive thresholding) and 
th e P X y (using FFT techniques) image registration estimators are 
the most attractive suboptimum techniques for small data image 
sizes (less than n=N=16). For large data image sizes, the p 

(using FFT techniques) has ’ computational advantages . 

c) In general, the radial error performance of the most widely used 
image registration estimator, p , was not as good as that obtained 
with the v=l and v=2 estimators. 

d) In general, image registration using hard limited (2 level) images 
provided greater radial error than that achieved using 64 level 
imagery . 

e) To determine the best parameters for data image size and shape or 
to determine which spectrum, terrain type, or image resolution is 
best for multispectral image registration will require a large 
number of runs to establish reliable statistical averages. 
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2.7 AREAS FOR FURTHER STUDY 

This initial study of digital image correlation techniques has suggested 
several interesting areas for further investigation. For example, there 
are two constraints on the optimum position location estimator which should 
be relaxed to make it more applicable to real world image registration pro- 
blems. These are: 1) the assumption that the difference image samples are 
statistically independent, and 2) the assumption of a uniform a priori 
probability of alignment. The removal of these constraints should be 
examined and their computational cost and effect on position location per- 
formance determined. This can be accomplished by extending the present 
optimum position location technique formulation. 

In addition, new position location techniques have been developed which 
have potential for becoming the "best" suboptimum technique. These recently 
developed techniques should be examined and included in the simulation for 
comparisons of computational efficiency and radial error performance. 

Finally, the most obvious area for further work is to perform a sufficient 
number of simulation runs to provide statistically reliable results. Image 
registration parameters such as data image size, shape, and quantization 
should be optimized for computational efficiency. Simulation results 
should then be obtained to determine radial error performance as functions 
of position location technique, multispectral image spectrum, terrain cate- 
gory, and image resolution. Fast search algorithms should be examined to 
determine their effectiveness at reducing computational cost. 
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3.1 INTRODUCTION 


SECTION 3 

IMAGE SAMPLING STUDY 


Previous studies have developed theories for image sampling and image recon- 
s true tion ^(References 10 and 11)* Little has been provided to date, however, 
which presents a quantitative comparison of techniques or which relates to 
actual image statistics. The purpose of this study is to provide a quantita- 
tive evaluation of several image sampling strategies and methods of image 
reconstruction. 


Two sample densities are considered in this study: (1) One sample for each 
58.06 square meters (625 square feet), and (2) One sample for each 3,716 
square meters (40,000 square feet). These densities bound the range of 
sampling densities encountered in practice* Image statistics based upon actual 
aerial imagery are utilized. These statistics were obtained in Section 1.4. 
Parametric studies based upon an average mean squared error of reconstruction 
are presented which show the effect of varying sample densities, terrain 
types, aperture sizes and shapes, and sample lattices. Also included is a 
measure of the errors introduced by presample filtering and image recording 
device characteristics. 

This study is confined to the accuracy of reconstruction for an image based 
on samples obtained through a sensing device which has an undistorted view of 
a scene. No image distortion correction is involved nor any time dependency 
considered. For purposes of the analysis, an infinite scene having specified 
spatial statistical properties is considered. 

3.2 IMAGE SAMPLING AND RECONSTRUCTION PARAMETERS 

Relevant parameters relating to image sampling and reconstruction are dis- 
cussed in the following sections. 

3.2.1 Sampling Lattice 

The sampling lattice is defined as the set of points (denoted by A) in the 
spatial domain (image plane denoted by X) which defines the centers of the 
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instantaneous fields of view of a sensor at the instant a particular sample 
is observed. Each of these points has a sample value of image intensity. 

This collection of samples is the data from which the image is to be recon- 
structed. Sample lattices in this study are usually sets of discrete points 
in two dimensions. They can form a square lattice, a rectangular lattice, or 
a rhombic lattice. These are defined in Figure 3-1, a, b, and c respectively. 
The spatial coordinates are (x^, x^) denoted by x* The names of the lattices 
refer to the shapes of a basic cell in the lattice. Even more general 
sampling strategies are possible as candidates but usually the strategy is 
determined by the shape of the sensor instantaneous field of view (IFOV) and 
these generally lead to the lattices mentioned above. 

Since continuously scanned images (video) are also of interest, we can include 
this possibility in the category of a sample lattice which is infinitely dense 
in the scan (x^) direction so that the set A in this case is a family of 
parallel lines . 

In connection with the sample lattice, there is another point set, B. This 
is a point set in the spatial frequency plane (ft) which is "dual" to the set A 
in the sense that if (a^ , a^) are the basis vectors (see Figure 3-1) for 
set A, then vectors (b^ , b^) in ^ satisfying 

a . .b . = 2 tt<5 , . (3-1) 

-i “j 

determine a set B in ft . This is the two dimensional set of folding frequencies 
introduced by the sampling procedure and is more fully explained in Reference 
10 and in Appendix A of this report. 

3.2.2 Sampling Aperture 

In general, a sensor will not be able to determine the image intensity at a 
single point, but will give some average value of intensities in the neighborhood 
of the point. This can be described by convolving the image with a weighting 
function y (x) corresponding to the sensor IFOV. For example, if the aperture 
is a uniform disk (that is, the sensing device senses every point in a circle 
with equal intensity and nothing outside the circle) , then the aperture 
function is as shown in Figure 3-2. Usually the volume under the surface 
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a. Square Lattice | 
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b. Rectangular Lattice 
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c. Rhombic Lattice I 


FIGURE 3-1 

IMAGE SAMPLING LATTICES I GP72-0450-3 


/VfCDO/V/VELL 4/f?CI7Arr 

3-3 \ 



REPORT MDC A1740 


defined by this function is taken as unity. 


%>| 



FIGURE 3-2 

CIRCULAR APERTURE WEIGHTING FUNCTION 


Other apertures of interest include square, rectangular, hexagonal and gaussian. 
The weighting functions for these are given in Appendix B, together with their 
two dimensional Fourier transforms, which will be of interest in the subsequent 
analysis . 

The size and shape of the aperture selected for best performance may depend on 
the sample lattice used, since too much overlap will smear the reconstructed 
image, and too small an aperture will omit information present in the image. 

3.2.3 Presampling Filtering 

Electronic or optical filtering can be performed on the signal representing the 
image before it is digitized and transmitted over the communication channel. 

In this study the effect of real time filtering of a video image is investi- 
gated. This is represented as a first order lag on the order of a fraction 
of a sample interval. 

3.2.4 Interpolation 

Image reconstruction by means of interpolation between sample values is 
studied by comparing three different interpolation methods. The simplest of 
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these and indeed the simplest of any interpolation scheme is one-point or step 
interpolation in which the reconstructed image at any point depends only 
on one sample value. Actually the reconstructed value is the value of the 
sample which is nearest in some sense. This is discussed in more explicit 
detail in Appendix C, as is the following method. 

An intermediate reconstruction method is a four point scheme which performs 
an interpolation on a sample cell determined by four surrounding points. 

This interpolation is a two dimensional generalization of common linear 
interpolation of one dimension. 

Also an interpolation is studied which is optimal in a least square error 
sense. This method for image reconstruction is theoretically the best that 
can be achieved with a particular aperture and sample lattice. This least 
squares estimate is best described as a two dimensional Wiener filter with 
no noise and no causality constraint and which depends on sampling aperture, 
sampling lattice and image statistics. It is derived in Appendix E. 

3.2.5 Recording Device Characteristics 

The output of the sampling and reconstruction procedure is not usually the 
(piecewise) continuous function determined by the interpolation method, but 
rather, a superposition of spots arrayed in a recording lattice. The shape 
of the spots depends on the recording device while their magnitude depends on 
the interpolated value at the spot center. The degradation in image reconstruc- 
tion introduced is examined in this study for two recording spot characteristics, 
a gaussian and a rectangular spot. The printing lattice is considered to be a 
refinement of the sampling lattice for analytical tractability in the study 
(that is, the sample lattice is a sublattice of the recording lattice). 


3.3 METHOD OF ANALYSIS 
3.3.1 Error Criterion 

The frequency or wave number space formulation of Peterson and Middleton 
(Reference 10) , which has been specialized to two dimensions, forms the basis 

i 

for most of the analysis in this study. The band limited condition has been 
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eliminated since the actual scene power spectral densities to which this study 
is directed do not exhibit this characteristic. An expression for the 
averaged mean square error of image reconstruction has been derived in 
Appendix A and is as follows (Equation A-31) . 


E{[f(x) - f(x)] 2 } 


( 2 ? 0 2 


$ (to) 


r (to) r (-to) 


I 1 ' i 
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u e B 
—m 


-u-uj j 


dw (3-2) 


where 


$ (tt>) = scene spatial power spectral density, 

T (u>) = Fourier transform of aperture function, 

G (oj) - Fourier transform of interpolation function, 

Q = area of sample lattice cell, 

Q, = two dimensional spatial frequency space, 

B = lattice in ft dual to sample lattice, 

= ~ spatial frequency (radians per meter) 


These quantities were more fully explained in Section 3.2. This expression 
is generalized to include filtering (Appendix G) and recording device charac- 
teristics (Appendix D) . It is this mean square error of image reconstruction 
which forms the primary basis for comparison of sampling strategies, apertures 
filtering, and recording characteristics in this study. The procedure for 
numerically evaluating these expressions is discussed in the next section. 


It should be kept in mind that averaged mean squared error is not the only 
possible basis which may be used for evaluation of sampling strategies. For 
particular applications (boundary detection, for example) one may wish to 
weigh certain spatial frequencies more than others. This can easily be 
accommodated by restricting the region of integration to a desired subset of 
ft or more generally by inserting a frequency dependent weighting function 
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into the integrand. Another technique which may be used is a two dimensional 
generalization of the modulation transfer function approach discussed in 
Appendix F. Results from these other approaches are, for the most part, 
beyond the scope of this study, but are not intrinsically difficult to in- 
corporate. 

3.3.2 Numerical Evaluation of Mean Square Error 

Evaluation of the expressions for mean squared error (Equation 3-2) must be 
done numerically. One-dimensional studies of sampling and reconstruction such 

i ■ -■ ' 

as the study of Ref erence 12 are much simpler to perform as integration over 

1 

the real line is required rather than over the entire plane. If in a one- 
dimensional case, the integrand functions are rational algebraic, then the 
theory of residues can be used for evaluation, greatly simplifying the 
procedure at no loss in accuracy. In two dimensions, however, this is not 
possible. The roll-off characteristics of the integrand functions are of 
primary concern in influencing the value of the mean squared error. This 
effect is much more pronounced in two dimensions, necessitating great care 
| in the numerical integration procedures. 

- Mean squared error computation requires numerical quadrature in which the 
^ integrand is a fairly complicated expression involving power spectral densities 
aperture and reconstruction functions, and lattice sums. Furthermore, the 
integration is over an infinite two-dimensional domain. Since the expressions 
involved are of a somewhat arbitrary nature, techniques such as gaussian 
quadrature are not applicable, and, since the expressions are of great com- 
plexity, high order methods involving partial derivatives are not practical. 
Because of this, the selected method is a trapezoided rule using Romberg 
quadrature (Reference 13) • Because the integrand is the difference between 

i ■ 

two comparable quantities, numerical integration error must be monitored. 

Romberg's method is a scheme by which higher order integration techniques 
can be generated by a linear combination of the trapezoidal sums with 
decreasing step size (h) . For the one-dimensional case, the trapezoidal 
integration T is an h^ process so that the value of the integration (y) is 
in error by an amount eh^. 
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y = T(h) + eh 2 


(3-3) 


If the integration step size is doubled, the error increases by 2 2 or 

y = T(2h) + e (2h) 2 (3-4) 

Using the correct combination of Equation 3-3 and 3-4, the error term can be 
subtracted out. 


4y = 4T(h) + 4 eh 2 

- y = T(2h) + e (2h) 2 (3-5) 

3y - 4T(h) - T(2h) 


y - T < 2h > ■- S(h) 

3 


(3-6) 


which is Simpsons rule. But Simpsons rule is an h 4 process; therefore: 


y = S(h) + eh 4 
y = S(2h) + e(2h) 4 


Continuing as before 

2 4 y = 2 4 S (h) + e (2h) 4 
y = S(2h) + e ( 2h ) 4 

(2 4 -l)y = 2 4 S (h) - S(2h) 

y = 2 - S - ( - h) - S(2h > = N (h) 
2 4 - 1 


(3-8) 


(3-9) 


This yields a Newton-Cotes technique, which is an h 6 process. Therefore, a 
general expression for the Romberg T-table can be developed. 


^ntj, (k+1) _ T (k) 

m-1 m-1 
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4 m - 
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( 3 - 10 ) 
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which gives the T-table 


T 

T 

T 

T 


0 

0 

1 

0 

2 

0 

3 

0 


L 1 

„2 


2 


where m = 0 
m = 1 
m = 2 
m 2 


Trapezoidal rule 
Simpson’s rule 
Newton-Cotes rule 
higher order Romberg methods 


(3-11) 


Looking at 
Ax 2 . 


the two-dimension cases, the h in Equation (3-3) is replaced by 


4 

y = T(Ax) + eAx 
y = T(2Ax) + e2 4 x 4 


S( x) 


2 T (Ax) - T(2Ax) 
2 4 - 1 


(3-12) 


Therefore, the general expression for the T function is simply modified by 
a power of two. 


T~ = 
m 


4 2m T (k+l) _ T (k) 
m-1 m-1 


, 2m 


- 1 


(3-13) 


The above scheme has been used extensively to evaluate the average root mean 
square error expression of Equation (3-2) . 


In the evaluation of the integrand, several lattice summations must be per- 
formed which are of the form: 


^2 f <!2 + 

u e* 

— m • 


(3-14) 
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where and are integers and b_^ and b^ are the dual basis vectors. The 
function is then summed by incrementing m^ and m^ in a sequence as shown in 
the Figure 3-3. 



This allows a finite number of lattice rings to be calculated which constitute 
a great savings in computer run time. The number of rings is dependent on 
the functions which are being summed over. 


The final output of the computer program is the average root mean square 
error normalized to the standard deviation of the power density spectrum.', 
and this output is printed in the T-table format of Romberg’s Quadrature 
method. 


We may capitalize upon the fact that lattice sums are involved in the inte- 
gration in order to transform the integration from an infinite to a finite 
domain without any appreciable increase in complexity. This is done simply 
by making the integrand a lattice sum and restricting the domain of the inte- 
gration to one sample cell. Further simplification is effected by taking 
notice of the evenness of certain functions involved. A functional block 
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diagram of the program is displayed in Figure 3-4. 
3.3.3 Image Power Spectral Densities 


Figure l-8c shows a power spectral density obtained from a discrete Fourier 
transform of the farmland scene of Fig. l-8a, digitized into a 128 x 128 
array of pixels with sample spacings corresponding to 30 meters. This result 
is plotted as a function of|w|in Figure 3-5, the data having been averaged 
circularly about the origin in the Q plane. Since the mean square error in 
image reconstruction depends directly on the spatial high frequency compon- 
ents, this curve must be extrapolated analytically to be useful for these 
studies. We therefore attempt to fit to this plot a density of the form: 


$(a>) 


4tto 2 D 2 (M-l) 
c 



(3-16) 


The roll off factor M is easily seen to be 3/2. The correlation distance 
for best fit was judged to be about 679 meters. This curve is also plotted 
in Figure 3-5. For a sample interval to be used as the standard of 
distance as is done in the computer simulation, the number of samples per 

correlation distance is 
D 


D ~ A 25 ft. 


679 meters/interval 


89 samples 


x 


.3048 meters correlation interval 


ft. 


and 


D = 


sample 

for the 1 sample/58.06 sq meter (625 sq ft) case 
679 


(3-17) 


= 11 samples/correiation interval 


200( . 3048) 

for the 1 sample/3716 sq meters (40,000 sq ft) case 


Urban area, mountains, and forest regions were also investigated. These had 
roll off factors of M = 1.02, 1.54, and 1.37 respectively. M cannot be less 
than unity or the integrals may not exist. This can be seen by constructing 
a two-dimensional Parsival relation to obtain the mean squared image intensity 
from its spectrum. This suggests that the urban scene may require a very high 
sample density for reconstruction or that higher order representation need De 
made for its spectral representations. These latter will require a more 
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FIGURE 3-4 

ROOT MEAN SQUARED RECONSTRUCTION ERROR FLOW DIAGRAM 
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extensive processing of real image statistics. 

Comparison of the above results does show a possibility for the classifica- 
tion of terrain types by their spectral characteristics. 

3.4 QUANTITATIVE RESULTS 

With computer implementation of the averaged mean squared error expression, 
sample densities and lattice shapes, sample aperture overlaps and shapes, and 
recording device characteristics were then evaluated. The comparison of 
alternate strategies is emphasized since absolute quantitative specification 
of rms error is not entirely meaningful except for very specific applications. 

The following sections discuss these effects. 

3.4.1 Influence of Image Spatial Frequency Content 

The average mean squared error obtained from any specific application depends 
strongly on the image power spectral density 0(aO. For band limited functions, 
it is possible to eliminate reconstruction error entirely if the sample density 
is high enough. It would follow, therefore, that image roll off M (see Section 
3.3.3) is an important parameter. Figure 3-6 illustrates the effects of varying 
M. The optimal reconstruction algorithm was selected in order to isolate the 
effect under consideration. Subsequent analysis was performed using the 
M = 3/2 value. Other values of roll-off are expected to shift the results in 
the manner indicated by the figure. Note that an overlap of about .6 will give 
near optimal results. This is about the overlap required to leave no area of 
the image unrepresented. Larger overlaps than .6 are seen to be non-detrimental. 

3.4.2 Sample Density 

Sample density variation was investigated using two basic sample densities 
of one sample per 58.06 square meters (625 sq ft) and one sample per 3716 
square meters (40,000 sq ft). Figures 3-7 and 3-8 show that overlap character- 
istics are little affected. A circular aperture with the rhombic sampling 
lattice was selected as being representative enough to display these effects 
since lattice shapes are not very influential. 
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Rhombic Sampling Lattice 
Uniform Circle Sampling Aperture 
Optimal Reconstruction 
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EFFECTS OF IMAGE SPATIAL FREQUENCY ROLL-OFF 
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FIGURE 3-8 

EFFECTS OF IMAGE SAMPLE SPACING 
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3.4.3 Sampling Lattice Effects 

Peterson and Middleton (Reference 10) have shown that, for band limited i 
functions, the densest lattices obtainable without overlap are optimal. Thus, 
the rhombic lattice is best for circular and hexagonal aperture cases, and 
either the square or rhombic lattice is equally good for the square apertures. 
These results are carried over to the non-band limited case by examining the 
circular aperture case. Figure 3-9 shows a comparison of the square lattice 
with the rhombic lattice. These lattices have equal spacing in the x ^> x 2 
directions. In the rhombic lattice, the samples in adjacent scans are out of 
phase, so to speak. As can be seen from the Figure, the effect is slight 
from mean squared error considerations. 
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3.4.4 Comparison of Reconstruction Methods and Aperture Shapes 

A comparison of the reconstruction methods discussed in Section 3.2.4 has 
been made for the apertures described in Section 3.2.2. This was done for 
the sample density of one sample per 58.06 square meters (625 square feet). 
Since the definition of sample overlap is different for each aperture shape, 
the various apertures were considered separately. Figures 3-10 through 3-13 
show these results. The differences between the circle, hexagon and square 
are seen to be minor. The gaussian aperture difference can also be judged 
as small if overlap for that case is modified to be twice the aperture stand- 
ard deviation divided by sampling interval. These results indicate that 
efforts to improve reconstructed image quality would be better spent on 
interpolation technique improvement rather than aperture shape optimization. 

3.4.5 Effects of Continuous Scan (Video) 


Continuous scanning and reconstruction of the continuous record can be repre- 
sented by allowing the sample density in the scan direction to become infinite 
This, as was mentioned before, is equivalent to a dual lattice becoming sparse 
Allowing the sample density in the scan direction to become infinite trans- 
forms the integration over a lattice cell into an integral over an infinite 
domain. To avoid program modifications and to also gain insight into this 
limiting procedure, cases with increasing scan direction density were investi- 
gated. 

Figure 3-14 shows the rms error for the case of a gaussian aperture with two 
different overlaps between scan lines. The .5 overlap represents the situa- 
tion in which the one sigma circles are contiguous between successive scans. 
The 0.14 overlapp represents the condition existing for the McDonnell Digital 
Image Processing System (DIPS) for which the distance between scan lines is 
1.5X10“3 inches and the one sigma radius of the scanning aperture^ is .21X10”^ 
inches . 
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INTERPOLATION TECHNIQUE COMPARISON FOR CIRCULAR SAMPLING APERTURE 



FIGURE 3-11 

INTERPOLATION TECHNIQUE COMPARISON FOR 
HEXAGONAL SAMPLING APERTURE 
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INTERPOLATION TECHNIQUE COMPARISON FOR 
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Ratio of Scan Direction Density to Cross Scan Density) 

FIGURE 3-14 

CONTINUOUS SCAN CHARACTERISTICS 


Apparently, an increase in sample density (along the scan direction) by a 
factor of seven is equivalent to continuous processing. An approximate 
twenty-five percent improvement results from processing video rather than 
equally spaced samples. 

Figure 3-15 illustrates the effect of varying the interscan distance as the 
scan direction sampling density increases while the two dimensional sample 
density is held constant. (i.e. 1 sample/625 sq. ft.). The sampling 

apertures were assumed to be Gaussian such that their two dimensional spot is 
a one sigma elipse whose axis are in a ratio corresponding to the ratio of the 
scanning densities. For the example shown, the elipses were assumed to 
be contiguous (i.e. the overlap is 1/2). This curve illustrates the image 
degradation caused by aliasing errors in the absence of any presampling 
filtering . 


AfCOO/V/VEL JL 4ll7CI74Pr 


3-22 




Normalized Average RMS Error - % 


REPORT MDC A1740 


Constant Sample Cell Area = 625 sq ft 
Gaussian Sampling Aperture with an 0.5 overlap 



Ratio of Scan Direction Density to Cross Scan Density | 

FIGURE 3-15f 

EFFECTS WITH CONSTANT SAMPLE CELL AREA 
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3.4.6 Image Recorder Effects 

An expression for mean squared reconstruction error including the recording 
device characteristics is derived in Appendix D. Two cases were investigated 
with this expression; 1) a Gaussian sampling aperture with a rhombic 
lattice and a Gaussian recording spot characteristic, and 2) a rectangular 
sampling aperture with a rectangular lattice and recording spot characteristic. 

The results for the first case, shown in Figure 3-16, indicate that the 
recording spot spacing need only be one— tenth of the sample spacing to 
effectively have a continuous recording situation. The rms error is 
approximately doubled by using no interpolation at all. 



FIGURE 3-16 

RECONSTRUCTION ERROR AS FUNCTION OF IMAGE 
RECORDING DEVICE CHARACTERISTICS 
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For the second case, the rectangular sampling lattice has the same spacing 
as the rhombic in one direction but the spacing normal to this is four times 
as large so that the sample density is 3/8 or .2165 times the density for the 
Gaussian case. The aperture was rectangular in shape with this same 4:1 
ratio of sides, was uniform, and of a size to make the sampled FOV contiguous. 
The recording characteristics had these same features. The error resulting 
from this situation is also shown in Figure 3-16,. Apparently only about five 
recording spots between every pair of sample points is sufficient to give 
the same performance as a continuous recording device. 

3.4.7 Presampling Filtering Effects 

A simple instructive example of filtering effects which may occur before 
sampling is a first order time delay acting on the video before sampling. 
Figure 3-17 illustrates the added degradation incurred as a result of the 
time delay times scan speed being various fractions 0 of the interscan 
distance. The situation depicted is the same as for the .5 overlap 
continuous scan examined in Section 3.4.5. This was analyzed according to the 



FIGURE 3-17 

EFFECT OF FIRST ORDER DELAY PREFILTER 
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methods described in Appendix G. First order reconstruction is displayed 
since the optimal reconstructor would remove nearly all of the effects of 
the filter. The mean square error is shown for those cases with and without 
an anticipated time shift in the recons true ted image. The amount of the 
shift used was equal to the spatial delay of the first order filter. This 
shift is not optimum, but is a reasonable value. To be optimum, the shift 
would be selected as a function of the scene spatial frequency spectrum, 
scanning aperture, and reconstruction method used. 

3.5 IMAGE SAMPLING STRATEGY SPECIFICATION 

The sampling strategy to be selected will depend upon the applications envis- 
aged for the reconstructed imagery. The results discussed so far are based 
on the unequivocable minimums rms error criterion. This is desirable for 
most purposes, and if a comparison of apertures and lattices is required, 
this is the most natural and general approach. However, special applications 
can take advantage of more specialized techniques. Some of the requirements 
on sampling and reconstruction imposed by special user applications are 
discussed in the following section. We begin however, by pointing out some 
of the conclusions obtained from the rms error criterion. 

3.5.1 Minimum RMS Error 

\ 

For the cases examined in this study, we may conclude the following: 1) The 

sampling lattice should be rhombic for the circular, gaussian or hexagonal 
sampling apertures. The lattice can be either square or rhombic for the 
square or rectangular aperture (the error difference is approximately three 
percent); 2) Varying aperture overlap within the range of .125 to 1.0 
results in fifteen percent variation in rms reconstruction error; 3) The 
most significant tradeoff parameter appears to be the interpolation reconr 
struction technique which, between the extremes of zero order reconstruction 
and optimal reconstruction, contributes thirty percent to the rms error. 

These considerations indicate that sampling strategy has less effect than 
reconstruction considerations, at least for cases in which the sampling pro- 
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cedure is unconstrained. The constrained situation (such as sample densities 
in the x^, x^ 9 being quite different) is expected to impose significant 
requirements on the aperture shape, but this evaluation has not yet been 
performed. One would select a rhombic sampling lattice if no mission 
requirements counteract this choice. Also, no requirement for other than 
a simple gaussian aperture was demonstrated unless the sampling procedure is 
quite anisotropic. 

3.5.2 Position Location 


Position Location is defined as the spatial registration of two images. 
These images may differ due to different sun angles, different FOV, and 
seasonal changes. Sampling and reconstruction of these images must be 
weighted toward features which are of comparable dimension to the required 
registration error bound. Figure 3-18 shows a proposed procedure for 
selecting sample densities and weighting functions. First, the image power 
spectral density and the noise impose a resolution limit which must be 



FIGURE 3-18] 

POSITION LOCATION CONSIDERATIONS 
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greater than 2tt / (registration requirement) if the registration accuracy is 
to be met. Secondly, the sample density must be set at a density high 
enough to reduce the effect of aliasing at this frequency to be within 
acceptable bounds (say .1 of image power at that frequency). Finally a 
weighting function is selected for image reconstruction which passes only 
those frequencies useful for registration. 

3.5.3 Feature Extraction 

Feature extraction is such a general term that only selected examples can be 
discussed. As far as sampling and reconstruction are concerned, only those 
aspects of feature extraction which pertain to image spatial properties are 
of interest. An important concept in feature extraction is multispectral 
species identification. This concept involves the recording of imagery in 
several spectral bands, and the investigation of the intensity relationship 
of the various electromagnetic frequencies which determines a point in a 
feature space (Ref .j 14) . Usually , high spatial resolution is not of great 
concern, but this is a choice to be made on the part of the user. An impor- 
tant consideration is that the intensity reconstruction error be of compa- 
rable or less magnitude than the natural variation expected within a species 
class. This is necessary to keep the point from varying beyond the limits 
of the appropriate species identification procedure. 

Boundary detection is another aspect of feature extraction. This can be 
investigated by restricting the region of the spatial frequency domain to 
the band (annulus in the two-dimensional case) occupied by the boundaries 
of interest. In addition, aperture shape may be used to advantage in 
boundary detection if the orientation of the boundaries is fixed relative 
to the scan direction. This particular effect could be analyzed in the rms 
error sense by means of an anisotropic image power spectral density which 
represents the orientation of the boundary lines. 
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3.5.4 Random Sampling 

Random sampling of band limited functions theoretically can be used to 
advantage for power spectral density determination. Beutler (Reference 15 ) 
shows that the power spectral densities of such functions are determined by 
samples taken randomly at an average sample density lower than the Nyquist 
density. The demonstration of this does not, however, provide a construc- 
tive method of determining the power spectral densities from the samples. 
Reconstruction of the actual signal rather than its statistics is a more 
difficult problem. The practicality of such signal recovery at less than 
the Nyquist density is questionable from "stability" considerations. Landau 
(Reference 16) defines stable sampling as sampling with the property that 
small changes in the sample values produce correspondingly small changes in 
the recovered signal. Thus, stability in this context is a type of con- 
tinuity property. Such a stability requirement is present in an earth 
reconnaissance program. 

Quantitative results for certain types of random sampling situations have 
been obtained. Leneman and Lewis (Reference 17) have exhibited mean squared 
reconstruction errors for certain types of one dimensional random sampling 
of signals * which need not be band limited. The sampling procedures con- 
sidered were Poisson sampling and periodic sampling with skips. The recon- 
struction methods were rather restrictive but some generalization may be 
possible to two dimensions and higher order interpolation methods. This 
approach could give results of some practical importance, however the syste- 
matic utilization of random sampling for channel bandwidth reductions is 
an open question. 
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3.6 CONCLUSIONS' 

The intent of this study was to obtain quantitative evaluation of some of 
the basic features of sampling and reconstruction rather than draw specific 
conclusions from these evaluations. We have, in terms of rms reconstruction 
errors , quantitatively described the effects of parametric variation of such 
things as sample density, terrain type, discreteness of recording device and 
presample filtering. It is inappropriate to draw specific conclusions from 
these, however in trade off areas such as sample lattice, sample aperture 
size and shape, we can say the following. The mean squared error criterion 
does not impose stringent requirements upon aperture shape or size as a 
function of sampling lattice unless the sampling is very anisotropic. Also, 
we may generally conclude that interpolation to a refinement of about five 
recording device spots (in each direction) per pair of sample values is 
sufficient. That is, no further refinement will improve the reconstruction. 

It is best to allow mission constraints or instrumentation considerations 
determine the sampling procedure and then use the criteria discussed in this 
report to evaluate the consequences of adopting or adapting that sampling 
method. The study has shown how these criteria can be expanded to include 
very specialized applications. 

3.7 AREAS FOR FURTHER STUDY 

This study has suggested several interesting areas which merit further 
investigation. For example, more anisotropic sampling lattice configura- 
tions (preferably those compatible with advanced mission or instrumentation 
constraints expected in earth surveys) should be examined. Also, adaptation 
of the rms error criterion to specialized tasks such as boundary detection, 
registration, or feature extraction may provide useful results. 
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Noise superimposed on the image, either before or after sampling, is a 
limiting factor in many image sensing applications • This can be easily in— 
eluded into the present rms error formulation. Optimal reconstruction 
formulae incorporating this superimposed noise does not appear to add any 
significant computational complexities. 

Another area needing further development is in the derivation of better 
reconstruction methods. This has been shown to be a promising area of in- 
vestigation which allows meaningful improvement in image reconstruction. 

For example, development of practical suboptimal reconstruction formulae 
based upon judicious truncation of the optimal formulae warrants considera- 
tion. A computer simulation of image reconstruction which can represent 
interpolation of up to the fourth order would be an invaluable extension to 
the present study. Certainly the implementation problem must be of 
immediate concern. If there are practical limitations to any of the recon- 
struction methods, it is best to have these identified as early as possible. 

The addition of real imagery to supplement any of the above mentioned areas 
of investigation is a desirable feature to include. The illustration of an 
image actually displaying a five or ten percent reconstruction error will 
provide a valuable addition to the study. 
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Appendix A 

DERIVATION OF MEAN SQUARED ERROR EXPRESSIONS 


These derivations essentially follow Peterson and Middleton (Reference 10) 
taking just the parts which are of interest in deriving the required 
expressions. The condition that the scene spectrum be band-limited has 
been omitted as not being general enough for real applications. 

Let A denote the set of sampling points in the image plane X. Thus if 
_a^, a _2 are the basis vectors for the sampling lattice", 

A = {n 1 _a 1 +n 2 a 2 : n^ f n 2 integers}. (A-l) 

Let B denote the set of points in the frequency plane Q which is dual 
to A in the sense that 

A 

B = +m 2— 2 : m i ,m 2 inte § ers ^> (A-2) 


and 


a. • b . = 2tt6 . . . (A-3) 
“1 ~3 ij 

Let f(x) denote the image to be sampled. Assume this image is adjusted 
so that the average value of the intensity is zero. (i.e., E{f(x)} = 0) and 
an auto- covariance function is given. 


E{f (x)f ) } “ K(z). 


(A- 4) 


Denote the sample aperture function by y ( x) • The £th sample is then 




k 


f (z)y(v 0 ~Z) d Z> 


V^A. 


(A-5) 
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If g(x) is^ the reconstruction (interpolation) function, the reconstructed 




image function (x) , is given by 


E 


f*(x) = ) g(x-V £ )iKv £ ) • 


(A- 6) 


v,EA 


The mean squared error in the reconstruction is 

E{ [f (x)-f * (x) ] 2 } - E{f 2 (x)-2f(x)f*(x)+f* 2 (x)} = K(0)-2 ^ ' 

v £ EA 


g(2£ - Z £ ) 


)</c K ^-y 


y)Y(v^-x) d y + 


E E 


vEA vGA 

— £ — m 


u. 


(A- 7) 


g (x-v £ ) g /v /y K (z _ z.) Y (z 0 ~Z) Y (v^-z) d^dz 


■H 1 — m 


In the single integral make the substitution z =X“Z> dzfd^j and in the 
double integral make the substitutions o_=^_- , jr_=_z-v^ d0=d£, di_=dz_. Then 


E{ [f (x)-f* (x) ] } = K(0)-2 


E 8( -- i> /' 


K (z) y (z- (x-v. ) ) dz_ 


V,SA 


E E g ^ )8( 2'V2jfV X ^(i-xJvC-^Y^-Vj-DdodT 


v t SA VJ=A 


= K(0) 


:0)-2 ^ ^ g(x-V £ ) ^K(_z)y(z-(x-v^ 


(A- 8) 


) )dz 


Z £ EA 


EE g (x-v^ ) S y^K(o-t_)Y(-a)Y(v k “T)d£dx. 


V^A V^A 
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Let a bar over the expression denote the expectation uniformly averaged 
over a sampling cell ^ . Let Q denote the area of so that 


Q = l®i 'xi 2 l’ 


(A- 9) 


then 


E{[f(x)-f *(*)]> = K(0) 


- ^ ^ (z.) Y (ST ) ) d z^dx 


v,GA 




J ^ ^ 'y ] g(2i-Z £ )g(2i-Z Jl _ Z k ) Jx ^K(a-T_)Y(-a)Y(v k -x)d£dxdx. 


(A-10) 


-£ EA ^k eA 


A summation over a sampling lattice repeats the integration over all the 
sample cells in the plane so we may drop the summation and integration over 
the entire plane. Thus 


E{[f(x)-f*(x)] 2 }= K(0) - | 


X 6 ® Jx 


K ( _z ) y (zf-x) dzdx 



g(x)g(x-v k ) 


VjGEA 

-4c 



K (o-_t) y (~o) Y (Vj“T_) dodxchc . 


(A-ll) 


Denote the lattice of delta functions corresponding to the sample points 
by 


A(x) 


s s<2r -* K 

vfA 


(A- 12) 
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The Fourier expansion of this is 



u6B 
— m 


(A- 13) 


The nth Fourier coefficient can be obtained in the usual way, 


f ix* u v r 

v "’ 4<2£)d2 '2j 


uGB 
— m 


u GB 
— m 


-ix*u 1X*U 

n m 

e e dx 


-X ‘.f 


-ix* (u -u ) 

— — n m , 

e dx 


1 .1 


= Z fo fo e 


-i(?a 1 +na 2 )( V u m ) 


Qd?dn 


u EB 
— m 


(A- 14) 


u GB 
— m 


dn 


Representing 


u , and u as 
— n — m 

u = n,b- + n 0 b 
— i n 1—1 2—2 

u = m-b- + m 0 b 0 . 
—m 1—1 2—2 


(A-15) 
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with n^, , n^, m 2 integers we have 


I. 


i . .1 

-iA (u -u )s 
e 1 d$ 


■ // 


-i27r(n 1 -m 1 )C 


d5 


*, n l “l 


' -i2n (n..-in ) 

,e -ll/l-i 2 Tr(n 


^ m^) ^ > n i^ m i 


= 6n^m^. 


1 

Similarly J ” ^ ^—2 


n 


dn = 6. 


n 2 m 2 


hence 


f -ix*u % » 

/ e n A(x)dx = 7 c Q 6 6 

•Aj, “ “ / ; ™ n 1 m 1 n 2 m 2 


£ 


u GB 
— m 


u GB 
— m 


c Q 6 =c Q 
m nm n 


(A- 16) 


(A-17) 



1 

Q* 


A(x) 



u fc 


(A-18) 


(A- 19) 
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Using this expression the averaged mean squared error becomes 

e{; 


[f (x)-f*(x) ] 2 } = K (0 ) - /gOO (z ) y ( z-x) dzdx 

J.E /. /. g (x) g (X~£_) A A K(a-r 


ip_-u (A-20) 

) Y (— c ) y (p_-T_) e dadTdxdQ_. 


uEB 
— m 


Let G(w) and T (w) denote the Fourier transforms of K(x) , g(x) and 

yOO respectively so that for the transform of a transformable function 


f^(x) we have 


y r -io)*2L 

X f 1 (x)e dx, 


(A-21) 


and also 


f 1 (x) = 


(2ir)‘ 


P l ( ^ )e d “- 


(A-22) 


We shall make use of a Parseval theorem which states that 


-i (0) = — /^(ujdu. 

(2ir) Jq 


(A-23) 


Now consider the function h(_£> defined by 


h(y.) = K(x) - j^g (x-y_) K(x)Y(z-x)dzdx 
— 2 ^ ^ Jyg (x~x) / x§ (x-£_) e J'yYiQri) J'yf' (-STl) y(.~o)da_djdpdx. 


(A-24) 


u EB 
— m 


Then 


E{[f(x)-f*(x)] 2 } = h(0). 


(A-25) 
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Rewriting the expression we obtain 

h(y) = K(y_) - ^ I ^K(z.)y(-(x-z.))dz 



g(-(y-x))dx 

y (p - x ) dk 


1P_.II 

g <£-£.) e 


.... .XI /XIJX Y( -2.> K tL-^ d 2 

u^B 
— m 

We now derive the Fourier transform of a convolution of the form 


(A-26) 

g(-(y-x))dx. 




1^*2. 


f 3 (x) = /x f l^ f 2^-"^ e d ^- 


F 3 (o)) 


/• J 

/ “103 X 

jo) = / X f 3 (x) e dx 


■ /x/x £ i 


(y)f 9 (x _ X) e e d^dx 


/* -i (oj-sO *x -iaji (x~x) 

Jy. Jy£ i(y) e f 2 (x-y)e dyd(x _ x) 


- A 


-i(uj-jz)*X /* -1 oj*£ 

(y)e dx / x f 2 (e_) e dp_ 


(A-27) 


(A-28) 


F 1 (oj-x)F 2 (^) • 


The Fourier transform H(u)) of h(x) is obtained by repeated application 
of the above formula. 


H(oj) = $(o)) - ^ 0(o))r(-o))G(-o)) 


+ 




(A-29) 


r(-(u-u ) ) ® (oj— u ) r (tn— u )G(o))G(-w) 

m m m — — 


uGB 
— m 
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A Parseval relation then furnishes the desired expression for the average 
mean squared error. 


E{[f(x)-f*(x)] } = 


7?f*[ 

E 


$ (oj) - (to) r (- 0 )) G (-uj) 


(A- 30) 


H r 7 $(a)-u )r(w-u )r(-(co-u ))G(w)G(-u)) 

pZ / j m — m m — — 

u6B 
—m 


da). 


Since the sampling lattice is infinite, this expression can be written as 


E{ [f (x)-f*(x)] } = 


^~2 L 

(2tt) J 


) I 1 - Q r (-W.)G (-oj) 


E 


+ -~T (to) T (— _co) y G (oo+u_)G (- (oH-u ) ) 

Q 

u £B 
— m 


da). 


(A- 31) 


MCDO/V/VELL 4/PCW4FT 


A- 8 . 



REPORT MDCA1740 


i 

Appendix B 
APERTURE TRANSFORMS 


Determination of the Fourier Transforms of the sampling apertures is a 
straightforward application of the necessary integral formulae. A few of the 
cases will be listed in this appendix for reference purposes. Each of the 
apertures y (x) has been normalized so that the volume under the surface they 
represent is unity. Some of the apertures are depicted pictorially in 
Figure B-l. 


Gaussian Aperture ; 



(B-l) 


It should be remembered that this is a two dimensional situation so that the 
volume inside a la cylinder is 39.3% of the total. 2 a and 3a represent 86.4% 
and 98.80% respectively. The transform of this function is also gaussian. 



Upaform Disk of Radius R : 


Y(x) 


r(o>) 



o, |k|>r. 

2J-| (|h>jR) 
|w|R 


(B-2) 


(B-3) 


(B-4) 


Here J^( ) is the Bessel function of the first kind of the first order and |to| 
is the vector magnitude of (i.e. |co| = v4o*co ). 
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Square with edge 2d: 


y(x) = 


~r—n> x ^ square, 
4<T 

0 , x 0 square, 


(B-5) 


r(co) = 


sin co^d sin u^d 
a)^d co^d 


and are the coordinates of to (i.e. to = (eo^,^)). 
Hexagon of altitude d : 


Y (x) = 


.mlji i , ,x £ hexagon, 
2/3? 

0 , x. hexagon . 


r(w) = 


1 1 


3 u+v 


, . /.os sin u . . /o . \ sin v 
sin (u+2v) — — h sm (2u+v) — — — 


where 


1 

3 

1, 


u 


sin 2v 


v 


/ . \ 2 " 

/ sm v \ 

V ’ / J’ 


sin 2u + | sin u 


2 cos v + 


sin v 


sin v 


, u^O and v^O, 
and u-fv^O, 

u=0 and v^O 


u^O and v=0, 

u+v=0, v^O 
u=v=0. 


v = — 


H% +a *h 


(B-6) 


(B-7) 


(B-8) 


(B-9) 
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FIGURE B-1 1 
SAMPLING APERTURES| 
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Appendix C 

RECONSTRUCTION TRANSFORMS 


The simplest reconstruction in one dimension is step type interpolation in 
which the value of the reconstructed function is equal to the value of the 
nearest sample. This "one-point" reconstruction has a simple generalization 
to two dimensional reconstruction. This is defined by specifying that the 
reconstructed image f* as a function of the samples ip be 

f*(x) - i K v £ ) (C-l) 

with v^ being the closest sample point in A to x in the metric determined by 
the sample lattice basis vectors a^, a^* More specifically if a point x is 
represented by 

x = aa^ + ga. 2 (C-2) 

then its "distance" from the origin is 

2 2 1 /? 

II x. .11 = ( a + 3 ) 7 (C-3) 

so that 

II2S.II = + (x -Jb 2 ) 2 ] 1 ^ 2 (C-4) 

in terms of the dual basis vectors Jb^, b^. 


It is easily verified that for this case the reconstruction function g(x) is 
defined by 


- g (x ) = 

in other words g(x) = 1 on a 
Figure C-l. 


|x«b^2.l <1r an< i 1 1 <7r 

0, otherwise, 

rhombic figure determined by a.^, 


(C-5) 
as shown in 
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FIGURE C-1 | 

ONE POINT RECONSTRUCTION GEOMETRY 




, t 

I 

■..A 



J 


\ 


GP72-045O-81 


In transforming the function g(x) we change coordinates from x^, x^ to a, 3. 
The Jacobian of this transformation is easily seen to be 


J = 


3x^ 3x^ 
3a 33 

9x 2 3x 2 
3d 33 


= Q = area of sample cell. 


(C-6) 


Thus the Fourier Transform of g(x) is 


G(to) = Q 






yl/2 

r 1/2 

- to.* (aa.^+3a^2> 

.1 .1 
sin -^o.a.. sin — <o.a„ 

/ 

/ 16 


dadg = Q — j— = — j— ^ - 

7 -1/2 

-1/2 


s 2^1* £2 


(C-7) 
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A two dimensional generalization of one dimensional linear interpolation 
is the bilinear form 


g(x) 


(l-|a|) (1- | 6 | ) , for | a | and | 3 | <1 
0, otherwise. 


(C-8) 


For this case g(x) depends upon the four nearest points, that is, the four 
points contained by the rhombic figure determined by the sample basis vectors 
as shown in Figure C-2. 



FIGURE C-2 

FOUR POINT RECONSTRUCTION GEOMETRY 


GP72-0450-78 ] 


In this case 



(C-9) 


These results bear an obvious similarity to the one dimensional case. 
Further generalization along these lines is possible. Also mixture between 
basis vector directions of different interpolation methods is easily done. 
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Appendix D 

DERIVATION OF MEAN SQUARED ERROR EXPRESSION WHICH INCLUDES 
THE RECORDER SPOT CHARACTERISTIC 


To include the recording characteristics into the formulation let f(x) be 
the recorded image function. Then 


f(x) = 


E 

vGC 

-q 


f*(v^)p(x-v q ), 


(D-l) 


where p(x) is the recorder characteristic; that the recorder, if printing 

a unit value at will produce an image p(x). C is the set of points at 
which the recorder prints. This set is assumed to form a periodic lattice. 
Inserting the expression for the reconstructed image f* as a function of the 
sample values we have 


f(x) = 


HE 


8( VV*(V 


vGC 

-q . 


V^A 


p(x-v ) 




v|EA -v^C 


z t eA v^c 


T(v t ) 


(D-2) 




Assuming that the lattice C is a refinement of the lattice A, (ACC), we 
can write 


f(x) = g(v r )p((x-v £ )-v r ,0 

/ — A ^ __ t — r> 


v £ CA vGC 




(D-3) 
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Define 


g*(x) 



v6C 


g(v r )p(x-v r ) . 


Writing this as 


a convolution we have 


g*(x) 


E 

v£C 

~r 

E 

u GD 



g (z) P (x~z ) 5 (z _ z r ) d Z 


1 

g(x)p(x-z)-e dy_, 


where S is the area of a cell of the printing lattice C, 
in the U plane with basis vectors reciprocal to those of 
transform of this expression we obtain 


G*(a>) 


1 

S 


G (or- Ug ) P ( w) 

u GD 


where G*(aj) and P(u>) are the transforms of g(x) and p ( 2 c) 
Noticing that 


f(x) = 


g*(x-v fc )iKv a ), 


whereas 


E 


f*(x) = 7 t g (x-v £ )^(v £ ), 
-l eA 
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(D-4) 


(D-5) 


and D is the lattice 
C. Taking the Fourier 

(D-6) 


respectively. 

(D-7) 


(D-8) 
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we find that the averaged mean square error in the total reconstruction 
including the recorder characteristics can be obtained by substituting g*(x.) 



(D-9) 


•(u+u,)) 

— Tin 


dto 


u eB 
“in , 


where 


G* (to) 


1 

S 



G(0Jr- U )P(0>) ’ 


u ED 


(D-10) 


For a realistic recording situation the overlap characteristics of ^the 
printing spots is not a parameter which may be varied but depends upon the 
spacing so that the average intensity in any area represents the scene being 
reproduced. To determine this relationship we consider the limiting case in 
which p(x) becomes narrow and the density of printer spots becomes large, 
Q* ( e., C becomes dense, D becomes sparse and S becomes small).* 

Recalling that the reconstructed image is 


f(x) 



f*(v q )p(x-Vq) 


vp 


(D-ll) 


then if 
so that 


f*(v^) = l for all v^GC we require that the average value of f(x) be 1 



where Y is a recorder cell 
and S is its area 


(D-12) 


/*fCOO/V/V^f_I_ AIPCPAPT 

D-3 



REPORT MDCA1740 


thus 



p (x-v^)dx=S. 


(D-13) 


This requirement sets the "gain" for the recorder* This can also be 


written as 



or as 


(D-14) 


P(0) = S. 


(D-13) 
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Appendix E 

OPTIMAL LINEAR INTERPOLATION FOR NON-BAND LIMITED FUNCTIONS 


For certain non-band limited functions an optimum linear interpolation can be 
identified similar to a Wiener filter except that for two dimensional appli- 
cations the causality requirement is not as meaningful so that spectral 
factorization is not involved. The approach is by a variational procedure 
as was Peterson and Middleton’s with the exception that their prefiltering 
considerations are not used. Consider the formula given in Appendix A for 
the mean squared error (not averaged) as a functional of g(x). 


E[g(x)] = K(0) - 2 £ g (x-vp jx K (*) Y ( z.- (x _ X £ ) ) d .z 


v £ e A 


(E-l) 


+ 


El 8 (x~ ) 8 (a - t ) y (-a) Y (v ^-r ) da_dx_. 


vjEA vGA 

— £ — K 


We simply require that for all admissible h (x ) 


lim 

e->0 


9E [g(x) + eh(x)] = 0. 

3e 


(E-2) 


This gives 


-2 


£/. 


X K (z ) y (z- (x - y ^ ) ) d£ h (x-v^ ) 


v,e A 


/j 8( ^rv 

v^A vEA 


/x /x K <2- 


(E-3) 


t^)y(-ct_)y(v -T)dgdTh(x-v ) 

K J6 


y ^ g(x-v a ) J y J'yY* (o~l) Y (-£) Y dadjrh (x-v^-v^ ) = 0. 


v^SA v^gA 
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Rearranging terms in summations this becomes 


-2 



K (z ) y ( z - ( x- v ^ ) ) d zh ( x- v £ ) 




V£k V=A 



K(£-t_)y (“£)y (v 


-t) dad th (x 


+ ^ ^ g«S-v 4 )-v,) 

v^A vjEA 




X K (£ _ X) Y (-£ ) Y (-v k ~£) d£dih ( 


Implying that we must have 


C O /* /* Y(v -r)+y (-V- 

/ X K (z ) y (z-x ) d£ = y ^ 8(2L"Z k ) y x y x K(£-O y(-£) 2 


vEA 


In terms of a summation over the B lattice, 

( z) Y (z-x) d z_ = ^ ^ ^ J^e ““g (x-zj) K (£-£) Y (“£ 


Y (z~_T )+Y (-£“1 


)- 


u EB 
— m 


The Fourier transformation of this is 


$ (to) r (-uj) 


= 1 V 1 

Q / > 2 


[$(ujrHi )+$(-co“‘U )]F((aH-u )F(-oi)“U )G(a>) 
m m — -m ~*m — 


u EB 
-m \ 


so that the optimal condition is given by 

$(u))F (-oj) 


G (a)) = 
opt — ' 


tE $ (-co-u ) ] r (w+u ) r (-co-u ) 

} m m — — m m 


Q ^ [H ~ 


u EB 
— -ra 


-V 

(E-4) 

x-v £ ) = 0. 


l) 

— d£d£. 
(E-5) 


.) 

— d£d£d£. 

(E-6) 


(E-7) 


(E-8) 
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provided that the denominator does not vanish which is rarely the case for non 
band limited functions and if this expression is inverse Fourier transformable. 

In any practical situation the scene power spectral density § (_w) is an 
even function so that 


$ (to) r (-to) 



ugB 
— m 


Inserting this into the expression for the average mean squared error 
(Appendix A) we have 


E{ [f (x)-f* (x) J } = 

— opt — 


— — , f t*(«: 

( 2 tt) JQ 


) - ^ $(w)r(-a))G opt (-a)) 


+ -t- / $(tdfu )r(wfu )r(-u)-u )G (cj)G (-u)]dco 

^2 m m' m opt — opt — — 


u GB 
-m 


t— / $(oj) [1 - ^r(-0))G (-U))]da) 

(2n) 2 Ja ~ Q - opt - - 


(E-10) 


- ^2 f 

(2tt) Jq 


4> ( co) 


$ (co) r (o)) r (-td) 


1 - 


E 

u£B 
— m 


$(orl-u)r (oorfu)r (— tO“U.) 


dto. 


This is the formula to be used in calculating the optimal reconstruction 
error. Note that for band limited scenes (support of $ (oj) is restricted to a 
finite region of the Q plane) and high enough sample density, then the sum- 
mation appearing in the expression will consist of evaluation at the point 

u = {0} only and the expression reduces to 0. 

— m — 
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Appendix F 

MODULATION TRANSFER FUNCTIONS 

The modulation transfer function concept in common usage for optical systems 
(Reference 18) can be adapted to spatial frequency applications in image 
^analysis. In this way the resolution of reconstructed images can be pre- 
dicted using communications engineering methods for introducing any of the 
image degradations or enhancements encountered in a space data retrieval 
situation* 

This procedure is based upon the transformation properties of power spectral 
densities. If $^(u>) is the (two dimensional) power spectral density of the 
input of a linear system with weighting function y (x) and r(oj) is the Fourier 
transform of y (x ) then the output power spectral density Oq(w) is given by 

* 0 (w) = r(w)r(-w)$.(w) (F-i) 

1/2 

It will be seen that | T (to) T (-to) | may be thought of as a modulation trans- 
fer function. It can represent optical systems, atmospheric attenuation or 
electrical filtering in a scanning situation (see Appendix 1 G) . 

A sampled and reconstructed image requires a further sophistication which is 
in common usage in communications engineering. If A is the set of sample 
points forming a periodic lattice (with cell area Q) , and if B is the fre- 
quency lattice dual to A, and if g(x) with transform G(w) is the interpolation 
formula as in Appendix A, then the reconstructed image power spectral density 
$*(to) is given by: 


$*(ay) 


ecu? s<-») y-» » (» + u > 

Q Q Z-I ~ r ” 


u£B 

-m 


(F-2) 


The summation represents the folding frequencies introduced by the sampling. 
One reservation must be made in introducing the power spectral density of the 
sampled and reconstructed image because of the following considerations. 

The local frequency content of the reconstructed image will depend upon the 
position within a sample cell. This is because the covariance, E[f*(x)f*(z) ] , 
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is not a function of _z-x only. If we consider this as a function of j 2 -x and 
z_, then the _z dependence is periodic with the sample lattice. If this _z 
dependence of the covariance function is averaged uniformly over a sample 
cell the resulting function of z-x can be transformed to give an "average" 
power spectral density. This is the procedure needed to give the above 
result for $*(co). 

The modulation M of a scene is defined in terms of the darkest intensity 
D and the lightest intensity L by 


M = 


L-D 

L+D 


(F-3) 


so that there is a simple relationship between modulation and power spectral 
density. The intensity I at a spatial frequency co is I = . 


The average intensity 2 (zero frequency) of the scene is 2 - ^T~> hence 



(F-4) 


whence 


M(o0 



(F-5) 


This is a two dimensional modulation figure. Resolution in photographic 
applications is frequently expressed in terms modulation by empirically 
setting the limit of the human eye at a modulation of .01. Thus resolution 
can be expressed in terms of spatial frequency bandwidth which may differ 
depending upon angular orientation due to anisotropic behavior introduced by 
the scanning or data processing. 


The optical considerations such as atmospheric attenuation and optical system 
modifications can be done in terms of modulation transfer functions modifying 
the modulation as usual or translated into communications theory terms involv- 
ing power spectral densities. 
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Appendix G 

ANALOG SCANNING AND FILTERING 


Analog scanning can be represented as a limiting case of a rectangular 
sampling lattice in which the sample density in one direction becomes infi- 
nite thus turning the sampling pattern from an array of points into a pattern 
of parallel lines. The scanning aperture giving coupling between scan lines 
can thus be taken into account by retaining the two dimensional formulation. 

The formulation is actually simpler in this case since the dual lattice B 
becomes sparse in one direction and in the limit is nothing but a linear 
string of isolated points representing the line scan frequency and its multi- 
ples. Thus fewer terms need to be retained in the lattice sums involved in 
mean squared error calculations. 


Electrical filtering of a scanned image acts upon only one scan line at a 
time and can be represented as an aperture function involving a delta func- 
tion for its cross scan spatial dependence. This function can then be cas- 
caded with the actual aperture function acting before the sampling (scanning) 
takes place. As an example consider a first order time lag. This has the 
weighting function 

_ *1 

1 1 *"* 0 ' f 

7 T e <$(x 9 ), x <0 

(G-l) 

0 , x ^>0 


where the distance 0 is the product of the time constant and the scan speed. 
The coordinate x^ is in the direction of a scan line. The Fourier transform 
of this is 


W(w) = 


1 

l+im^0 


(G-2) 


which is complex but will only appear in the formulae as 


Re {W(o0 } = W(m)W(-m) = 2 ~ 1 ~ ’ 

l+u 1 0 


(6-3) 
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The mean squared error formula as modified to accommodate this is 



— m 


(G-4) 
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Appendix H \ 
NEW TECHNOLOGY 


There were no reportable items during the performance of this contract. 
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